Plot Wavelet FFT in Hz

Hi! I'm trying to do fft of wavelet function. But I'm confused: how to define x-axis f in Hz? Here is my code:

scales=1:100;   
wname='db20';   
N=10;            

[phi, psi, tgrid] = wavefun (wname, N);

f=???

figure;
hold on
plot(f,abs(fft(phi)));
plot(f,abs(fft(psi)));

before, in my practice, i used:

f=0:Fs/(N-1):Fs;

What should I do now?

Answer 1

[Robert U]

Hello Alexander,

if you want to display wavelet analysis in terms of frequency one common way is to use the center frequency of the used wavelet.

Using a web search engine you could have searched for "wavelet equivalent frequency" which is providing several useful explanations as: https://de.mathworks.com/videos/understanding-wavelets-part-1-what-are-wavelets-121279.html

Kind regards,

Robert

[Alexander Voznesensky]

Thank you. I tried:

[freq,xval,recfreq] = centfrq(wname, N, 'plot');

And I know the points number:

N=length(xval);

But what about Fs? How to define f vector?

Answer 2

[Robert U]

Hi Alexander,

if you want to plot the Fourier spectrum of the wavelet you can apply the usual fft algorithm (described here: https://de.mathworks.com/help/matlab/math/basic-spectral-analysis.html ).

 %%Create wavelet
 wname='db20';   
 N=12;            

 [phi, psi, tgrid] = wavefun (wname, N);

 %%Calculate FFT of wavelet output
 SampleFrq = 2^N;        % calculate sample frequency of signals
 SigL = length(psi);     % phi & psi have same length

 NFFT = 2^nextpow2(SigL);    % standard fft calculation
 SigMagPsi = fft(psi,NFFT)/SigL;
 SigMagPsi = 2*abs(SigMagPsi(1:NFFT/2+1));

 SigMagPhi = fft(phi,NFFT)/SigL;
 SigMagPhi = 2*abs(SigMagPhi(1:NFFT/2+1));

 f = SampleFrq / 2 * linspace(0,1,NFFT/2+1);

 %%Plot Figure
 figure;
 hold on
 plot(f,SigMagPhi);
 plot(f,SigMagPsi);
 ah = gca;
 ah.XScale = 'log';
 ah.YLabel.String = 'Magnitude';
 ah.XLabel.String = 'Frequency [Hz]';
 legend(['Phi';'Psi'])

 figure;
 centfrq(wname,N,'plot')

[Alexander Voznesensky]

Thank you very much! But, the main question is still here: why

SampleFrq = 2^N?

I have some stuff about it:

Fs=10^3;                  % частота дискретизации, Гц
T=1/Fs;                   % период дискретизации, сек
N=1e4;                    % количество отсчетов
t=(0:N-1)*T;              % массив отсчетов времени
f=0:Fs/(N-1):Fs;          % массив частот

f1=10;
f2=25;
f3=50;
f4=100;

x1=cos(2*pi*f1.*t)+cos(2*pi*f2.*t)+cos(2*pi*f3.*t)+cos(2*pi*f4.*t); % стационарный сигнал

figure;
subplot(2,1,1);
plot(t,x1);
xlabel('t, сек');
ylabel('x');
title('Стационарный сигнал. Временная область');
grid;

subplot(2,1,2);
plot(f,abs(fft(x1)));
xlabel('f, Гц');
ylabel('X');
title('Стационарный сигнал. Частотная область');
grid;

x2_1=cos(2*pi*f4.*t(1:N/4));
x2_2=cos(2*pi*f3.*t(N/4+1:N/2));
x2_3=cos(2*pi*f2.*t(N/2+1:3*N/4));
x2_4=cos(2*pi*f1.*t(3*N/4+1:end));

x2=[x2_1 x2_2 x2_3 x2_4]; % нестационарный сигнал

figure;
subplot(2,1,1);
plot(t,x2);
xlabel('t, сек');
ylabel('x');
title('Нестационарный сигнал. Временная область');
grid;

subplot(2,1,2);
plot(f,abs(fft(x2)));
xlabel('f, Гц');
ylabel('X');
title('Нестационарный сигнал. Частотная область');
grid;

So my Fs is 1000 Hz. What in this case? Is Fs the same for the Fourier spectrum of the wavelet? So:

f=0:Fs/(length(psi)-1):Fs;

[Robert U]

Hello Alexander:

The basic idea is to understand the generated wavelet as a time-dependent signal that is sampled with a certain (fixed) sampling frequency even though it could be a function depending on any physical (scalar) value (e.g. distance).

Sampling frequency is the step size of your points grid XVAL (tgrid). Assuming that they are equidistantly distributed you can check the equivalent sampling frequency by

 Fs = 1 / mean( diff( tgrid ));

You will see that the sampling frequency increases with the number of iterations. You can calculate your sampling frequency either with the formula above or you convince yourself that the scaling is according to

 Fs = 2^N;

This feature was documented on earlier versions ( http://radio.feld.cvut.cz/matlab/toolbox/wavelet/wavefun.html ) and is still for wavefun2 ( https://de.mathworks.com/help/wavelet/ref/wavefun2.html).

[PHI,PSI,XVAL] = wavefun('wname',ITER) returns the scaling and wavelet functions on the 2^ITER points grid XVAL.

Kind regards,

Robert