How to obtain Fourier power spectrum of velocity data?

Hello,

I would like to know which frequencies are hidden in my velocity signal (=velocity vector).

I started straight forward: I took the "fft" of the vector and plotted its amplitude against the frequencies, but this does not give me the desired result, i.e. no clear frequency peaks were observed.

Since the velocity vector is non-periodic, i tried to smoothen the beginning and end of the signal by means of a hanning window. Unfortunately, it did not work.

I attached all the data in the zip-file it contains:

• UpV_irreg.mat (velocity vector)
• t_irreg.mat (corresponding time vector)
• fft_test.m (my code)

Any help would be greatly appreciated.

Answer 1

[Star Strider]

See if the pwelch (link) function will do what you want. It may also not show any specific frequency peaks (it tends to combine closely-spaced peaks), however it will do essentially everything you describe as what you want to do.

[Lukos]

Thanks for the quick reply.

However, with the pwelch() command a similar (undesired) result is obtained.

When I use the Curve fitter tool and use a sinus function fit on the signal, I was able to obtain a frequency, but this is only a rough estimate.

So I figured, since the curve fitter tool is able to give me some kind of fit, that a fourier spectrum analyses would give me a more accurate frequency.

[Star Strider]

My pleasure.

The fft function is the most appropriate function to understand the spectral components of your signal.

I get reasonably good results with this code:

Dt = load('t_irreg.mat');
Dd = load('UpV_irreg.mat');
t = Dt.t_irreg;
D = Dd.UpV_irreg;

L = numel(t);
t = linspace(min(t), max(t), L);
D = resample(D,t);
Ts = mean(diff(t));                                                         % Sampling Interval
Fs = 1/Ts;                                                                  % Sampling Frequency
Fn = Fs/2;                                                                  % Nyquist Frequency
FT_D = fft(D-mean(D))/L;
Fv = linspace(0, 1, fix(L/2)+1)*Fn;
Iv = 1:numel(Fv);
[Pks,Loc] = findpeaks(abs(FT_D(Iv))*2, 'MinPeakHeight',3E-3);               % Fourier Transform Peaks & Locations (Index)
expstr = @(x) [x.*10.^floor(1-log10(abs(x)))  floor(log10(abs(x)))];        % Used For The ‘sprintfc’ Or ‘compose’ Call
EPks = expstr(Pks);                                                         % Used In The ‘sprintfc’ Or ‘compose’ Call
EFrq = expstr(Fv(Loc)');                                                    % Used In The ‘sprintfc’ Or ‘compose’ Call

figure
plot(t, D)
grid
xlabel('Time')
ylabel('Amplitude')

figure
plot(Fv, abs(FT_D(Iv))*2)
grid
xlim([0 0.1])
txtc = sprintfc('\\leftarrowAmpl = %4.2f\\cdot10^{%2d}\n    Freq = %4.2f\\cdot10^{%2d}', [EPks'; EFrq']')   % Choose One
txtc = compose('\\leftarrowAmpl = %4.2f\\cdot10^{%2d}\n    Freq = %4.2f\\cdot10^{%2d}', [EPks'; EFrq']')    % Choose One
text(Fv(Loc), Pks, txtc, 'FontSize',7, 'HorizontalAlignment','left', 'VerticalAlignment','top')
xlabel('Frequency')
ylabel('Amplitude')

I use the linspace and the resample functions to provide more regularly-spaced sampling (sometimes this can be important), then subtracted the constant (d-c) offset from your data before doing the fft. The plot presents the important peaks in your data, and annotates the amplitudes and frequencies at which they occur. The amplitudes are given in the ‘Pks’ variable, and the frequencies are ‘Fv(Loc)’.

[Lukos]

Thanks for helping me out! One more question regarding the frequencies I obtain now.

I realized my code gives me the same results, but you made me realize that the range of frequencies is in between 0 and 0.1.

Could you help me with one more thing:

The highest amplitude is at a frequentie f=0.03. Does this correspond to T = 1/f = 33.33?

Again, your help is much appreciated!

[Star Strider]

My pleasure.

‘The highest amplitude is at a frequentie f=0.03. Does this correspond to T = 1/f = 33.33?’

Not really. The highest amplitude (0.0183) is actually at 0.0278, and 1/0.278 = 36. (I did these calculations with the full-precision results.)

Calculate the periods of the peak frequencies as:

T = 1./Fv(Loc)'  

If my Answer helped you solve your problem, please Accept it!

[Lukos]

Thanks!

[Star Strider]

As always, my pleasure!

[Lukos]

How do I interpret the amplitude in the frequency domain?

I would like to convert the amplitude in the frequency domain to the time domain, i.e. what is the amplitude of the velocity data (in the time domain).

The amplitudes in the frequency domain where obtained by normalizing the fft of the veloicty signal with the number of samples L and by multiplying its magnitude by a factor 2 (since it is one-sided?). What scaling do I need to perform to obtain the time domain amplitudes of the velocity signal? Thanks in advance!

[Star Strider]

The amplitude in the frequency domain and time domain are the same in the code I posted. Note that MATLAB computes a two-sided Fourier transform, so the energy is equally divided between the positive and negative frequencies. It is therefore necessary to multiply the amplitude by 2 in the plotted one-sided Fourier transform to display the amplitude accurately:

plot(Fv, abs(FT_D(Iv))*2)

The discrete nature of the fft usually results in some ‘leakage’ of signal energy to near-by frequencies, so the peak amplitude in the fft may not always match (may be less than) the amplitude of the time-domain signal.

If you want to calculate the time-domain signal from the frequency-domain signal, you need to use the original two-sided Fourier transform results, without the amplitude scaling necessary for the one-sided Fourier transform.

I do not understand your reason for doing this, however. If you want to do frequency-selective filtering of your signal, do that in the time domain with some of the filter designs that the Signal Processing Toolbox (among others) allows you to create.