8

I am using Cmor-Fb-FC wavelet transform in matlab environment. Fb is bandwidth and Fc is center frequency parameters. with trail and error procedure I select Fb-Fc as 5-1, and I've got rational output. But I want to know whats the rule of selecting proper value for Fb-Fc parameters? Imagine I am using 300 samples over 30 sec. so my sampling frequency is 10 Hz.Here is my signal.

x=sin(2*pi*t*.1).*(t<10)+sin(2*pi*t*0.3).*... 
    (t<30)+sin(2*pi*t*0.6).*(t<10).*exp(-t*.1);
SAH
  • 1,317
  • 4
  • 17
  • 39

1 Answers1

10

Remember that Wavelet Transforms are nothing but time-localized filtering/correlation operations. The wavelet transforms provide a unified framework for getting around the Heisenberg Uncertainly Principle that the Fourier Transform suffers from. So when you ask "what should my settings be for bandwidth, and center frequency", you are asking for filter parameters to be given to you. No one can give you filter parameters but yourself, and your application. This is why they are parameters. :-)

More concretely, how to pick your $F_c$ and $F_b$ have to do with the nature of your signal. First, some pre-liminaries:

Wavelet Transforms are nothing but matched filters that you construct, such that you have the best match to the 'features' in your signal of interest. One set of 'features' that you might consider would be the neighbourhood of center frequency, and the neighbourhood of bandwidth. Again, no one can pick those for you, because you need to know what you are looking for. However, you can pick a range commensurate with what you expect. (This is actually one of features that make Wavelet Transforms very powerful).

The morlet wavelet is given by the following function:

$$ \Psi(t,f_c,T_p) = \frac{1}{\sqrt{\pi T_p}} \ e^{\frac{-t^2}{T_p}} \ e^{-j2 \pi f_ct} $$

Here, $t$ is the time index in seconds, $f_c$ is your center frequency, and $T_p$ represents what I have called a "Period Parameter". Inspect the above equation. Notice how the morlet wavelet is in fact nothing but a complex exponential centered at frequency $f_c$, windowed by a zero-mean gaussian function, with $\sigma = \sqrt{\frac{T_p}{2}}$

Thus, to answer your first question, the $F_c$ parameter, must be chosen such that it is at or near the frequency you wish to interrogate.

Regarding the wavelet bandwidth: What I have called $T_p$ here is your "Period Parameter", and this is directly related to your bandwidth. Go back to the Morlet Wavelet equation, and consider a scenario where instead of a gaussian window, we used a rectangular window. In the Fourier Domain, we attain a sinc function due to this type of window, and we can say that the Null-To-Null Bandwidth, $B$, is simply: $B = \frac{2}{T_{period}}$, (in the passband).

The figure below illustrates this for a complex exponential located at $f_c = 2$ Hz, with varying $T_{period}$ values:

enter image description here

You can clearly see that as the period extent is increased, the bandwidth decreases. This is the inverse nature of the time-frequency relationship. If you wanted to use this pseudo-wavelet, what would you pick the $T_{period}$ to be? The answer would be "Whatever $T_{period}$ my signal-of-interest hangs around at". In fact, this is exactly what is done in the Short-Time-Fourier-Transform, (STFT), although one would pick a smarter window instead of the boxcar.

Now let us re-inspect the Morlet Wavelet Equation, and re-insert the gaussian windowing function, with $\sigma = \sqrt{\frac{T_p}{2}}$. If we do this, we get the following time-domain plots, and the corresponding Discrete Fourier Transform (DFT)s of the Morlet Wavelets:

enter image description here

Once again, notice how as the extent - as the variance of the gaussian window shading the complex exponential increases, the bandwidth in the Fourier Domain decreases. (In this case, we have much nicer side-lobe behavior, but this is tangential). In any case, You can once again see how we have an inverse relationship between the $T_p$ parameter, (or the standard deviation parameter of your gaussian, whichever suits you), and the bandwidth of the corresponding filter.

What to do with this knowledge: Now you understand how bandwidth is related to the gaussian shading function embedded within the Morlet Wavelet. Thus, for any signal $x[n]$ you have, you may simply take the DFT of your signal, and observe its bandwidths at various center frequencies. Then, adjust the $\sigma$ and $f_c$ parameters of your Morlet Wavelet accordingly, such that the DFT function of the Morlet filter directly overlaps the DFT function of your template signal. This is how you can set your parameters.

The code to generate the above plots is here:

clear all;
fS = 500;
tStart = -4;
tStop = 4;
timeVector = linspace(tStart,tStop, (tStop-tStart)*fS );
fCenter = 2;
tP = [0.2 1 3];
exps = length(tP);

figure(1); clf(1);
figure(2); clf(2);
for ii = 1:exps

    tPeriod = tP(ii);
    timeMask = zeros(1,length(timeVector));
    timeMask((timeVector >= -tPeriod/2) & (timeVector <= tPeriod/2)) = 1;

    psiFakeWavelet = exp(2*1i*pi*fCenter.*timeVector).*timeMask;
    psiWavelet = ((pi*tP(ii))^(-0.5)).*exp(2*1i*pi*fCenter.*timeVector).*exp(-timeVector.^2/tP(ii));

    %Demonstrating the Rectangular Window
    figure(1);
    subplot(3,2,2*ii-1);
    plot(timeVector,real(psiFakeWavelet), timeVector,imag(psiFakeWavelet));
    ylim([-1.2 1.2])
    xlabel('Time / Seconds');
    title(sprintf('Morlet Wavelet, T_{period} = %2.2f, f_c = %2.2f Hz', tPeriod, fCenter));

    input = psiFakeWavelet;
    Nfft = 10 * 2^nextpow2(length(input));
    psd = 20.*log10(fftshift(abs(fft(input,Nfft))));
    freqs = [0:Nfft - 1].*(fS/Nfft);
    freqs(freqs >= fS/2) = freqs(freqs >= fS/2) - fS;
    freqs = fftshift(freqs);
    figure(1);
    subplot(exps,2,2*ii);
    plot(freqs, psd); 
    xlim([-10  15]);
    xlabel('Frequency / Hz');
    title (sprintf('PSD, Null-Null-BW = %2.2f Hz', 2/tPeriod));

    %Demonstrating the Morlet Wavelet
    figure(2);
    subplot(3,2,2*ii-1);
    plot(timeVector,real(psiWavelet), timeVector,imag(psiWavelet));
    ylim([-1.2 1.2])
    xlabel('Time / Seconds');
    title(sprintf('Morlet Wavelet, T_{period} = %2.2f, f_c = %2.2f Hz', tPeriod, fCenter));

    input = psiWavelet;
    Nfft = 10 * 2^nextpow2(length(input));
    psd = 20.*log10(fftshift(abs(fft(input,Nfft))));
    freqs = [0:Nfft - 1].*(fS/Nfft);
    freqs(freqs >= fS/2) = freqs(freqs >= fS/2) - fS;
    freqs = fftshift(freqs);
    figure(2);
    subplot(exps,2,2*ii);
    plot(freqs, psd); 
    xlim([-10  15]);
    xlabel('Frequency / Hz');
    title ('Wavelet PSD');

end
Tarin Ziyaee
  • 2,831
  • 16
  • 14
  • thanks! so if i have 3 modes in my signal([0.1,0.3,0.6]Hz ), then my Fc will be 0.25 Hz (0.6-0.1)/2 ? and how can I determine bandwidth of my signal?(Fb) – SAH Aug 05 '13 at 16:42
  • "getting around the Heisenberg Uncertainly Principle" I don't think that's quite correct. They still have the same problem, but they make different trade-offs at different frequencies. – endolith Aug 05 '13 at 19:10
  • 3
    @endolith The Fourier Transform contains 100% frequency information, and 0% time information, courtesy of the HUP. Wavelets get around this via a trade-off. – Tarin Ziyaee Aug 05 '13 at 19:12
  • @endolith,@User4619. thanks for information, would you answer my question about selecting the fb and fc please? – SAH Aug 05 '13 at 19:18
  • @Electricman I am editing the post. – Tarin Ziyaee Aug 05 '13 at 19:19
  • @user4619: Oh right, I'm thinking of STFT vs wavelets, not FT vs wavelets – endolith Aug 05 '13 at 20:24
  • 1
    @Electricman: It depends on the application. Try different values and see what highlights the features you're interested in. – endolith Aug 05 '13 at 20:37
  • @endolith, I did try. and I've got nice result for my test signal[the function I wrote on the post.]. but I was looking for a formula to determine Fb-Fc.cause my actual signals are not that abstract function, they are power system oscillations. And thanks for your comments. :) – SAH Aug 05 '13 at 20:53
  • @endolith Yes, that is exactly correct. – Tarin Ziyaee Aug 05 '13 at 21:04
  • @user4619 hi,Very useful information.I've tried to plot those wavelet and their DFT (or PSD) as you did. but I've faced to a problem with getting DFT. I will post my code and plot in the answer below. Please look and tell me whats the problem with my code at the DFT bit. thanks in advance – SAH Aug 06 '13 at 11:42
  • @Electricman Please see my edit with the code. – Tarin Ziyaee Aug 06 '13 at 13:34
  • @user4619 Very helpful dude. kind regards. – SAH Aug 06 '13 at 18:38
  • @user4619 - how is the frequency domain version of the morlet wavelet that you exposed? – user7780 Feb 07 '14 at 14:37
  • @user7780 What do you mean? I do not understand your q. – Tarin Ziyaee Feb 07 '14 at 14:51
  • @Tarin Ziyaee This is great material. Any chance you could apply to the OP's original signal? There are three frequencies of interest. Would you put the central frequency at one of them and scale to multiply/divide and get the other two? Would bandwidth then be chosen in order to resolve these three frequencies? Knowing we have to resolve these modes without leaking, how do we get the most time resolution possible? – Eli S Aug 27 '21 at 17:38