How to scale the frequency axis after performing fft?

Question

Hasan Alkhater 2016-9-14

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/303075-how-to-scale-the-frequency-axis-after-performing-fft

评论： ASAD RASHEED 2019-1-16

Hello, I am trying to analyse signals aquired via an accelerometer for a person walking. suppose I choose the sampling frequency to be 100 Hz. then I would scale the frequency vector as follows

 % code
N = length(Xabs);
fgrid = fs*(0:(N-1))/(N);

After I plot, the x-axis of the plot is scaled based on the sampling frequency being 100 Hz. Say the location of the dominant frequency in the plot is 4Hz. Now, if I change the sampling frequency to 1000, the location of the dominant frequency is ten times the previous location.

I have a feeling that my explanation is kinda messy. Here is the question: How can I scale the frequency axis such that, no matter what the sampling frequency is, the location of the dominant frequency is correct and does not change?

here is the complete script I use for your reference

fs = 100;
X = fft(x); % Obtain the DFT using the FFT algorithm
Xabs = abs(X); % Obtain the magnitude
N = length(Xabs);
fgrid = fs*(0:(N-1))/(N);
Xabs = Xabs(1:floor(N/4));
fgrid = fgrid(1:floor(N/4));
plot(fgrid,Xabs);

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

ASAD RASHEED 2019-1-16

Why are you taking the faxis and magnitude till N/4? Should it be till N/2?

请先登录，再进行评论。

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

Answer 1

dpb 2016-9-15

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/303075-how-to-scale-the-frequency-axis-after-performing-fft#answer_234876

在 MATLAB Online 中打开

Plot versus frequency instead of samples...the dominant frequency won't necessarily be in the same place but it'll show up at the correct frequency (and the same input signal sampled at differing rates will be the same location, of course)...

L=length(y);                      % series length
Fs=100;                           % or whatever is actual sampling frequency
f = Fs/2*linspace(0,1,NFFT/2+1);  % single-sided positive frequency
X = fft(y)/L;                     % normalized fft
PSD=2*abs(X(1:L/2+1)))            % one-sided amplitude spectrum

If you apply the accelerometer sensitivity to your input signal, you should get actual acceleration with due respect for no windowing, etc., ...

8 个评论
显示 6更早的评论隐藏 6更早的评论

dpb 2016-9-15

编辑：dpb 2016-9-15

在 MATLAB Online 中打开

You didn't do something correctly, then. Show your work...

Example:

>> Fs=100;  % sample @ 100 Hz
>> T=10;    % collect data long enough for at least a couple cycles
>> N=T*Fs;  % compute number samples needed for above
>> t=linspace(0,10,N);  % and generate time vector
>> y = 0.7*sin(2*pi*4*t)+randn(size(t));  % and a sample signal around 4Hz
>> Y = fft(y)/N;              % FFT
>> PSD=2*abs(Y(1:L/2+1));     % and the PSD one-sided
>> f=linspace(0,Fs/2,N/2+1);  % compute freq vector for Fs
>> plot(f,PSD)                % plot the result
>> hold all                   % we're going to put another on top...
>> [~,ix]=max(PSD)  % and where's the peak located???
ix =
  41
>> f(ix)
ans =
   4
>> Fs=250;  % ok, let's change sample rate
>> N=T*Fs;  % next steps just repeat for new sampling vector...
>> t=linspace(0,10,N);
>> y = 0.7*sin(2*pi*4*t)+randn(size(t));
>> Y = fft(y)/N;
>> PSD=2*abs(Y(1:N/2+1));
>> f=linspace(0,Fs/2,N/2+1);
>> plot(f,PSD)
>> xlim([0 20])  % move scale down to show peak area more clearly
>>

The resulting plot is

Note here I modified the sample size to ensure generate long enough time series to get at least a couple cycles in--this makes the two df the same for the two cases but as can be seen, they're different lengths but the full-width Fmax isn't the same.

You've got to ensure you're using the right sample times and such to get the normalization right. The key relationships are

Fmax=1/2dt (Nyquist)
T = N dt   (Total time)
df= Fmax/(N/2) --> 1/T (frequency resolution)
T = 1/df

dpb 2016-9-17

在 MATLAB Online 中打开

>> x=load('xacc.mat');
>> x
x = 
  xacc: [1401x1 double]
>> x=x.xacc;
>> t=load('time.mat')
t = 
  time: [1x1401 double]
>> t=t.time;
>> t(1)
ans =
   0
>> t(end)
ans =
   7
>> t(2)
ans =
  0.0050
>> 1/t(2)
ans =
 200
>>

AHA! There's the error--sampling rate is 200 Hz if this is indeed 7 sec of time history as dt shows...if it were 100 Hz sampling rate the total T would be

>> length(x)/100
ans =
 14.0100
>>

So, since it must actually be 200 Hz sampling rate, then

>> dt=1/200;
>> N=length(x);
>> Fmax=1/2/dt  % it's the Nyquist that's 100 Hz, NOT sample rate.
Fmax =
 100
>> df=Fmax/(N/2)  % so the frequency resolution is
df =
  0.1428
>> Fs=200;
>> f=Fs/2*linspace(0,1,N/2+1);
>> X=fft(x)/N; PSD=2*abs(X(1:fix(N/2)+1));
>> plot(f,PSD)
>> xlim([0 10])
>>

The resulting PSD has peaks at a number of frequencies between 2 and 8 Hz...

And, to show it works out correctly,

>> f(2)
ans =
  0.1429
>>

so the computed df is the same as that shown to be required above from Fmax. Again, if you apply the above relations to the actual sampling rate and overall time duration you'll end up at the correct frequency for the result; the relationships are the same as that I used in the demo for a given known frequency; you simply apply whichever are those that are the givens for a particular data acquisition to compute the other. Here it is assumed the total time is really 7 sec and so the dt for the number of samples enforces what the sample rate must have been. From that total time we also get what the max frequency can be as 1/T and that, with the sample length, produces the frequency bin resolution. IOW, to get higher resolution in the frequency domain, you must sample longer time for baseband analysis; otherwise must use zoom or other analysis techniques.

dpb 2016-9-18

编辑：dpb 2016-9-19

在 MATLAB Online 中打开

Glad to try, anyways... :) There's a lot to pick up, agreed; the above relationships of sampling are key and will (eventually) become second-nature.

I'll make one additional note on a practical matter in implementation/computation -- there's an "off-by-one" error easily made; for example the theoretical of

T = N dt

for the time history length; if you use a Matlab vector t which is always from 1:N, length(t) is N (of course). But there are only N-1 dt increments so the actual end time is (length(t)-1)*dt. That's where the 0.01 extra seconds came from in the case I used in the example with your data files to illustrate that a 100 Hz sample rate would be a 14-sec time history instead of the expected 7 as was 14.01 instead of 14 exactly. You have to be careful in computing the df for the same reason; note that there was a little difference in that precise numeric value between the approximate one of 0.1428 I showed as what the value needed to be as opposed to the end result of 0.1429 from the actual computed frequency vector.

That difference is the difference between an exact 1/7 which rounds off to the 9 in the fourth digit as opposed to

>> Fmax/(N/2)
ans =
    0.1428
>> Fmax/((N-1)/2)
ans =
    0.1429
>>

It doesn't make much difference on an individual value of df, but that difference will accumulate over a number of channels and eventually can lead to difficulties in analyzing frequency peaks if it's need to compare harmonics of bearings or the like, say, by causing enough shift that the ratios don't show up as they should.

Hasan Alkhater 2016-9-20

Got it I was thinking that this is because it is an estimate, which is why there is this slight difference

but it turns out to be due to the actual number of increments not being used

thank you ^^

dpb 2016-9-22

Yeah, I was too lazy to type the extra parens around the (N-1) term initially and the approximation was good-enough to prove the point of where the error lay...

请先登录，再进行评论。

How to scale the frequency axis after performing fft?

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

采纳的回答

8 个评论
显示 6更早的评论隐藏 6更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

How to scale the frequency axis after performing fft?

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

采纳的回答

8 个评论 显示 6更早的评论隐藏 6更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

8 个评论
显示 6更早的评论隐藏 6更早的评论