A couple of Time Frequency Analysis MATLAB programs

Mark DeArman

Mark DeArman

Junior
Aug 15, 2011
389
0
16
Puyallup, Washington, United States
www.cascadeacoustic.com
These should work in MATLAB or Octave.

A simple pseudo Wigner-Ville Distribution.
Code: 
[FONT=courier new]function [ WD, f, t ] = wvd( x, Fs )[/FONT]
[FONT=courier new]%WVD compute discrete wigner-ville distribution[/FONT]
[FONT=courier new]%   x - analytic signal[/FONT]
[FONT=courier new]%   Fs - sampling frequency[/FONT]
[FONT=courier new]
[/FONT]
[FONT=courier new][N, xcol] = size(x);[/FONT]
[FONT=courier new]if N<xcol[/FONT]
[FONT=courier new]    x = x.';[/FONT]
[FONT=courier new]    N = xcol;[/FONT]
[FONT=courier new]end[/FONT]
[FONT=courier new]
[/FONT]
[FONT=courier new]WD = zeros(N,N);[/FONT]
[FONT=courier new]t = (1:N)/Fs;[/FONT]
[FONT=courier new]f = (1:N)*Fs/(2*N);[/FONT]
[FONT=courier new]
[/FONT]
[FONT=courier new]for ti = 1:N[/FONT]
[FONT=courier new]    taumax = min([ti-1,N-ti,round(N/2)-1]);[/FONT]
[FONT=courier new]    tau = -taumax:taumax;[/FONT]
[FONT=courier new]    WD(tau-tau(1)+1,ti) = x(ti+tau).*conj(x(ti-tau));[/FONT]
[FONT=courier new]end[/FONT]
[FONT=courier new]
[/FONT]
[FONT=courier new]WD = 2.*fft(WD);[/FONT]
[FONT=courier new]
[/FONT]
[FONT=courier new]end[/FONT]

An inverse FFT which always returns real values.
Code: 
function [ x ] = idft_real( X )
%Return inverse FFT always real
Y = real(X) + imag(X);
y = ifft(Y);
x = (real(y)+imag(y));
end

An alternate to the MATLAB hilbert function.
Code: 
function [ z ] = analytic( s )
%ANALYTIC Compute the analytic function


N = length(s);
S = fft(s,2*N);
S(2:N) = S(2:N).*2;
S(N+2:2*N)=0;
z = ifft(S);
z(N+1:2*N) = 0;


end
 
Re: A couple of Time Frequency Analysis MATLAB programs

Sorry I guess I should have put the example of how to use it here. All of these functions need to be placed in files which have the name of the function and then all dropped into the same directory in octave.
There are two more functions which are useful for making log sweeps.

Generate the phase of a log sweep, this is from the integral of the ideal IF.
Code: 
function [ tt ] = logphase( t,T,w1,w2 )
%LOGPHASE Logrithmic phase trajectory


tt = ((T*w1)/log(w2/w1))*(exp( (t/T)*log(w2/w1)) - 1);


end

Generate the log sine sweep
Code: 
function [ x t ] = logsweep( f1, f2, T, A, Fs )
%LOGSWEEP generate sweep from w1 to w2 of length T
%   f1 - sweep start frequency in Hz
%   f2 - sweep end frequency in Hz
%   T - duration of sweep in seconds
%   Fs - sampling rate in samples per second
%   A - amplitude




t = 0:1/Fs:(T-(1/Fs));


w1 = 2*pi*f1;
w2 = 2*pi*f2;


tt = logphase( t,T,w1,w2 );


x = A.*sin(tt);


end
% end sinesweep

Example of how to use the WVD to visualize the sweep.
Code: 
function example
%EXAMPLE 
f1 = 20;
f2 = 200;
T = 1.0;
A = 1.0;
Fs = 1000;
x = logsweep( f1, f2, T, A, Fs );
s = analytic(x);


[ S, f, t ] = wvd( s, Fs );
contourf(f,t,abs(S));
end
 
Re: A couple of Time Frequency Analysis MATLAB programs

Martyn ferrit Rowe said:
Hey Bennett,
I found this useful :-
http://www.amazon.com/Processing-Wavelets-Electrical-Engineering-Computer/dp/0977858200

ferrit
Click to expand...

I would like to point out that the WVD and other Time Frequency Representations are not from the field "Wavelets" for the most part. Although the two fields of study are closely related. In fact the way wavelets are used in the Wavelet TFR is not typical, they are used more as a window function. Not that there are not lots of cool things that wavelet processing can do. But you will turn up a lot more useful information on Time Frequency visualization and processing if you don't search for wavelets.
 
Top Bottom