matlab FFT. Stuck understanding relationship between frequency and result

Question

We're trying to analyse flow around circular cylinder and we have a set of Cp values that we got from wind tunnel experiment. Initially, we started off with a sample frequency of 20 Hz and tried to find the frequency of vortex shedding using FFT in matlab. We got a frequency of around 7 Hz. Next, we did the same experiment, but the only thing we changed was the sampling frequency- from 20 Hz to 200 Hz. We got the frequency of the vortex shedding to be around 70 Hz (this is where the peak is located in the graph). The graph doesn't change regardless of the Cp data that we enter. The only time the peak differs is when we change the sample frequency. It seems like the increase in the frequency of vortex shedding is proportional to the sample frequency and this doesn't seem to make sense at all. Any help regarding establishing a relation between sample frequency and vortex shedding frequency would be greatly appreaciated.

Answer 1

Take a look at this related question. While it was originally asked about asked about VB the responses are generically about FFTs

Answer 2

this is probably not a programming problem, it sounds like an experiment-measurement problem

i think the sampling frequency has to be at least twice the rate of the oscillation frequency, otherwise you get artifacts; this might explain the difference. Note that the ratio of the FFT frequency to the sampling frequency is 0.35 in both cases. Can you repeat the experiment with higher sampling rates? I'm thinking that if this is a narrow cylinder in a strong wind, it may be vibrating/oscillating faster than the sampling rate can detect..

i hope this helps - there's a 97.6% probability that i don't know what i'm talking about ;-)

Answer 3

The problem you are seeing is related to "data aliasing" due to limitations of the FFT being able to detect frequencies higher than the Nyquist Frequency (half-the sampling frequency).

With data aliasing, a peak in real frequency will be centered around (real frequency modulo Nyquist frequency). In your 20 Hz sampling (assuming 70 Hz is the true frequency, that results in zero frequency which means you're not seeing the real information. One thing that can help you with this is to use FFT "windowing".

Another problem that you may be experiencing is related to noisy data generation via single-FFT measurement. It's better to take lots of data, use windowing with overlap, and make sure you have at least 5 FFTs which you average to find your result. As Steven Lowe mentioned, you should also sample at faster rates if possible. I would recommend sampling at the fastest rate your instruments can sample.

Lastly, I would recommend that you read some excerpts from Numerical Recipes in C (<-- link):

• Section 12.0 -- Introduction to FFT
• Section 12.1 (Discusses data aliasing)
• Section 13.4 (Discusses FFT windowing)

You don't need to read the C source code -- just the explanations. Numerical Recipes for C has excellent condensed information on the subject.

If you have any more questions, leave them in the comments. I'll try to do my best in answering them.

Good luck!

Answer 4

If it's not an aliasing problem, it sounds like you could be plotting the frequency response on a normalised frequency scale, which will change with sample frequency. Here's an example of a reasonably good way to plot a frequency response of a signal in Matlab:

Fs = 100;
Tmax = 10;
time = 0:1/Fs:Tmax; 
omega = 2*pi*10; % 10 Hz
signal = 10*sin(omega*time) + rand(1,Tmax*Fs+1);

Nfft = 2^8;
[Pxx,freq] = pwelch(signal,Nfft,[],[],Fs)
plot(freq,Pxx)

Note that the sample frequency must be explicitly passed to the pwelch command in order to output the “real” frequency data. Otherwise, when you change the sample frequency the bin where the resonance occurs will seem to shift, which is similar to the problem you describe.

Answer 5

Hi again. I tried using the frequency response code as above but it seems that I dont have the appropriate toolbox in Matlab. Is there any way to do the same thing without using fft command? So far, this is what I have:

   % FFT Algorithm

Fs = 200;                     % Sampling frequency
T = 1/Fs;                     % Sample time
L = 65536;                    % Length of signal
t = (0:L-1)*T;                % Time vector
y = data1;                    % Your CP values go in this vector

NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2);

% Plot single-sided amplitude spectrum.
loglog(f,2*abs(Y(1:NFFT/2))) 
title(' y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')

I think there might be something wrong with the code I am using. I'm not sure what though.

Answer 6

A colleague of mine has written some nice GPL-licenced functions for spectral analysis: http://www.mecheng.adelaide.edu.au/~pvl/octave/

(Update: this code is now part of one of the Octave modules:
http://octave.svn.sourceforge.net/viewvc/octave/trunk/octave-forge/main/signal/inst/.
But it might be tricky to extract just the pieces you need from there.)

They're written for both Matlab and Octave and serve mostly as a drop-in replacement for the analogous functions in the Signal Processing Toolbox. (So the code above should still work fine.)

It may help with your data analysis; better than rolling your own with fft and the like.

Answer 7

Methinks you need to do some serious reading on digital signal processing before you can even begin to understand all the nuances of the DFT (FFT). If I was you, I'd get grounded in it first with this great book:

Discrete-Time Signal Processing

If you want more of a mathematical treatment that will really expand your abilities,

Fourier Analysis by Körner