%Program wvd.m

%signal
r=8000; %sampling rate
t=(0:2047)/r;
omega1=2.0*pi*500;
omega2=2.0*pi*750;
f=exp(i*omega1*t.*t)+exp(i*omega2*t.*t.*t);

t1=0.192; %location of the delta function
t2=0.196;

k1=t1*r;
f(k1)=f(k1)+6;
k2=t2*r;
f(k2)=f(k2)+6;
plot(t,f)
figure(2)

%wvd computation
N=length(f);

if(mod(N,2) ==1);
   f=[f 0];
   N=N+1;
end

N2m1=2*N-1;
Nb2=N/2;

for m=1:N
    s=zeros(1,N2m1);
    s(N-(m-1):N2m1-(m-1))=f;
    s=conj(fliplr(s)).*s;
    s=s(Nb2:N2m1-Nb2);
    shat=abs(fft(s));

%Normalize with the number of overlapping terms
 
    if m<=Nb2
        shat=shat/(2*m-1);
    else
        shat=shat/(2*N-2*m+1);
    end
    
 wvd(m,:)=shat(1:Nb2);
end

%Plot

time=(0:N-1)/r;
freq=(0:Nb2-1)*r/N/2;

NColor=256;
colormap(gray(NColor));
brighten(-0.99)
wvd_min=min(min(wvd));
wvd_max=max(max(wvd));
wvd=(wvd-wvd_max)*NColor/(wvd_min-wvd_max);

image(time,freq,wvd');
%Because of the finite support of the signal, there will end effects

xlabel('Time (seconds)');
ylabel('Frequency (Hz)');
set(gca,'YDir','normal')
clear all