FFT frequency range/sample frequency simple problem

36 次查看（过去 30 天）
Michael
Michael 2013-10-11
评论： Michael 2013-10-12
I'm writing a code to process vibration acceleration time-history data by filtering into 1/3-octave bands, and I need to check that the filtered data is being processed correctly. To do this I would like to compare with a general FFT of the the time-history data. However, when I've tried to code this, I've run into problems getting the Matlab FFT algorithm to do what I want. So I've reverted to processing simple sine waves in an attempt to get a correct answer from that, before going back to the more complex acceleration data.
I've isolated the FFT part of the code, which is based on the example code found in the fft function information. However, my understanding of the finer points of shifting the data and setting the transform length and frequency range is still quite poor.
This is the isolated part of the code as it currently stands:
t = linspace(0,1,1000);
fs = length(t)/(max(t) - min(t));
y = 0.5*sin(2*pi*100*t) + 0.5*sin(2*pi*50*t);
m = length(t); % Window length
tl = pow2(nextpow2(m)); % Transform length
P1fft = fft(y,tl); % DFT
%P2fft = fft(P2,tl); % DFT
f = (0:tl-1)*(fs/tl); % Frequency range
P1ffts = fftshift(P1fft); % Zero-shift periodogram
%P2ffts = fftshift(P2fft); % Zero-shift periodogram
% FFT
figure
% semilogx(f,abs(P2ffts),f,abs(P1ffts))
plot(f(1:floor(tl/2)),abs(P1ffts(1:floor(tl/2))))
grid
% xlim([min(f) max(f)])
% xlim([2.5 100])
% ylim([0 0.014])
% set(gca,'XTick',Fref)
I note that if I set the sample frequency to 300Hz, this gives the correct output, but any other sample frequency does not. Clearly this must be related to the frequency range, but my brain is not able to make what is no doubt a fairly simple connection between these two bits of information, which is due to my lack of understanding of the how the code is running.
I have searched through the FFT threads already existing for some help, but perhaps this problem is so simple most people don't need to ask!
Thanks in advance for any help,
Mike
  2 个评论
Michael
Michael 2013-10-12
I have investigated this further and found that the issue is in relation to the fftshift function; taking this out (instead using only fft), and rejigging the code resolves the problem.
I had thought that fftshift was a necessary step in displaying a meaningful fft but this is clearly not the case.
The functioning code is now:
t = linspace(0,1,1000);
fs = length(t)/(max(t) - min(t));
N = (max(t) - min(t))*fs;
y = 0.5*sin(2*pi*100*t) + 0.5*sin(2*pi*50*t);
m = N; % Window length
nfft = pow2(nextpow2(m)); % Transform length
f = (0:nfft/2-1)*(fs/nfft); % Frequency range
P1fft = fft(y,nfft); % DFT
P1fft = P1fft(1:nfft/2); % Discard unnecessary part of spectrum
figure
plot(f,abs(P1fft))
I am still having some problems getting the amplitude scaling correct though - any tips would be well-received!
Cheers,
Mike
Michael
Michael 2013-10-12
PS. Any clear advice on correct usage of the fftshift function would also be appreciated - the help documentation is fairly opaque to me on this...

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Michael
Michael 2013-10-12
I have investigated this further and found that the issue is in relation to the fftshift function; taking this out (instead using only fft), and rejigging the code resolves the problem.
I had thought that fftshift was a necessary step in displaying a meaningful fft but this is clearly not the case.
The functioning code is now:
t = linspace(0,1,1000);
fs = length(t)/(max(t) - min(t));
N = (max(t) - min(t))*fs;
y = 0.5*sin(2*pi*100*t) + 0.5*sin(2*pi*50*t);
m = N; % Window length
nfft = pow2(nextpow2(m)); % Transform length
f = (0:nfft/2-1)*(fs/nfft); % Frequency range
P1fft = fft(y,nfft); % DFT
P1fft = P1fft(1:nfft/2); % Discard unnecessary part of spectrum
figure
plot(f,abs(P1fft))
I am still having some problems getting the amplitude scaling correct though - any tips would be well-received!
Cheers,
Mike
  1 个评论
Michael
Michael 2013-10-12
PS. Any clear advice on correct usage of the fftshift function would also be appreciated - the help documentation is fairly opaque to me on this...

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Fast Fourier Transforms 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by