Azzera filtri
Azzera filtri

Wigner Ville distribution - Integrating in time domain

4 visualizzazioni (ultimi 30 giorni)
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')
  3 Commenti
Sveinung Attestog
Sveinung Attestog il 21 Set 2020
%% 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
Sveinung Attestog il 21 Set 2020
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

Accedi per commentare.

Risposte (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by