Wigner Ville distribution - Integrating in time domain

Have been trying to make by own function for computing the Wigner Ville distribution. One consern I have is that my estimation seems to have simular shape, but the amplitude is smaller by a large factor (In this case a factor of 2000).

n=1000;
endtime=1;
time=linspace(0,endtime,n*endtime);
Noise=normrnd(0,0.1,[1,endtime*n]);
xk1 = sin(2*pi*50*(time));
xk2 = 0*sin(2*pi*49*(time));
xk=xk1+xk2;


%% My Wigner Ville distribution algorythm
%Hilbert transform of data set
hil_xk=hilbert(xk);
%Defining time parameter
tau_delta=endtime/(n);                      %Increments of tau during integration
t=endtime*0.6;                              %Set time which I'm computing for 
t_start=0;                                  %Start of data set
t_end=endtime;                              %End of data set
tau_max=2*min(t-t_start,t_end-t);           %Computing the limits for tau which is necesary for integrating. Either x or x* will be 0 if it's any larger
%Setting up frequency which will be computed
f=(0:1/(2*endtime):100-1/(2*endtime));
tau=(-tau_max:tau_delta:tau_max)';
%Computing all element of the function that will be integrated
x=interp1(time,hil_xk,t+tau/2);
xconj=conj(interp1(time,hil_xk,t-tau/2));
Xj=xconj.*x.*exp(-2*pi*1i*tau*f);
%Integration with resect to tau for different values of frequency f
Q=trapz(tau,Xj)/(2*pi);


%The signal
figure
plot(time,xk)

W2=wvd(xk');

f2=linspace(0,n/2,(n*endtime));
figure;
plot(f2,W2(:,endtime*n)*0.1/200)
xlim([0 100])
hold on
plot(f,(real(Q)'));
legend('Matlab','My function')
grid
title('Wigner distribution t=0.5')

[Manisha Dixit]

Hi..How will it work if we need to calulate WVD for an image(which is a matrix of order "lets say- 300*300,),without using the inbuilt function wvd.

[Sveinung Attestog]

%% Numerical example
fs=1000;
T_start=0;
T=5;
Time=linspace(T_start,T,fs*T);

Xk1 = sin(2*pi*25*Time);
Xk2 = sin(2*pi*60*Time);
Xk3 = sin(2*pi*20*Time.^2);
Xk=Xk1+Xk2+Xk3;

%% Genrating Time-Frequency distribution

sigma=1e-3;
tic
Cplot=Cohen(Xk,T,fs,sigma);
toc
%% Plotting the results
Ntetta=linspace(0,(fs)/2,fs*T);
time=linspace(0,T,T*fs);

figure
imagesc(time,Ntetta,Cplot)
ylim([0 300])
colorbar
set(gca,'YDir','normal')
xlabel('Time (s)')
ylabel('Frequency (Hz)')
%caxis([0 0.01])
%xlim([0 40])
colormap((parula))

function Cplot=Cohen(xk,T,fs,sigma)
hil_xk=hilbert(xk);

%% Auto Correlation function
Ntau=((1:1:fs*T)-T*fs/2-1);
Nt=(((1:1:fs*T)')-T*fs/2-1);
hil_xk=[hil_xk 0];

n1=flip(Nt)+Ntau;
n2=flip(Nt)-Ntau;

%Filtering out indexes outside of sample array

n1(n1>fs*T/2-1)=T*fs/2;
n2(n2>fs*T/2-1)=T*fs/2;
n1(n1<-fs*T/2)=T*fs/2;
n2(n2<-fs*T/2)=T*fs/2;

n1=n1+fs*T/2+1;
n2=n2+fs*T/2+1;


x=hil_xk(n1);

xconj=conj(hil_xk(n2));

Ax=(x.*xconj);

%% Ambiguity Function

A=fftshift((ifft(Ax)),1);

%% Kernal Function
%Chaning eta and tau 0-value to close to 0 value for avoiding NAN+NANi 
Neta=linspace(-fs/2,fs/2-1,fs*T)';
Ntau=Ntau*2/fs;
Neta(Neta==0)=1e-9;
Ntau(Ntau==0)=1e-9;

%Choi Williams Kernal function changed to a Wigner Ville kernal function
alpha=1/sigma;
P_CW=1;%exp(-alpha*(Neta*Ntau).^2);

%% Charateristic function
M=P_CW.*A;

%% Infers Fourier over eta
% t=0 results that the exponetial factor of the function is always 1 
%, thus we can execute the invers Fourier just by summation
C1=fft(fftshift(M,1));

%% Fourier transform
%Aquiring the final frequency dist. 
C=fft(C1');

Cplot=flip(((flip(abs(C)))'))';
end

[Sveinung Attestog]

Hope this will help. It's one of my scipt, which I probably based on another script I found online. It was currently set to a Choi-Williams kernal function, but I changed it to Wigner-Ville for you. Hope this will help. This example generates a matrix 5000*5000 for the images (Time period: 5 sec, sampling frequency: 1 kHz). It comutes in less than 2 sec (Checked by matlabs tic toc function).

Sorry if it is a little bit messy