How to scale the frequency axis after performing fft?

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);

[ASAD RASHEED]

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

Answer 1

[dpb]

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., ...

[Hasan Alkhater]

I think that the NFFT was meant to be L?

Anyways, same problem occurs. Once I plot PSD vs f, the location of PSD values relative to frequency will change if I change the sampling frequency Fs for example, Fs is initially 100, the location of the most energetic frequency is around 4 Hz. Once I change Fs to 1000, the location changes to 40Hz. That means I can't rely on this plot to determine the frequency, or the bandwidth of the energetic frequencies.

If I try to estimate the dominant frequency from the time domain, it is around 0.9 Hz. Which was not true for both values of Fs after plotting the magnitude spectrum.

[dpb]

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

[Hasan Alkhater]

Okay.. I am going to explain what I am doing because I think it is different from generating a sine wave based on what the sampling frequency is. suppose you have xacc data and created a suitable time series, see the attachment. (the x component of the acceleration). xacc contains 7 seconds of the x component of the accelerometer reading. How would you process and plot the frequency domain of this data such that you have the right frequencies in place regardless what the sampling frequency is? again, here is the script I run (first I make x=xacc):

 fs = 100;            %the sampling frequency
N = length(x);        
X = fft(x)/N;            %Performing fft
PSD = 2*abs(X(1:N/2+1)); %Creating the PSD vector
fgrid = fs/2*linspace(0,1,N/2+1); %frequency vector
plot(fgrid,PSD);        %plot PSD vs frequency in hertz

[dpb]

>> 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.

[Hasan Alkhater]

okayyyy. looks like there is that much that I don't know about signal processing

Thank you for all your help!

[dpb]

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]

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]

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...