help-octave
[Top][All Lists]

Re: Any FFT expert around?

 From: Brian Kaczynski Subject: Re: Any FFT expert around? Date: Thu, 13 Mar 2014 17:51:03 +0100

Hi Oxy,

Not quite.  See below.

2014-03-13 14:57 GMT+01:00 oxy :
Hi Brian,

> freq=linspace(-Fs1*(N-1)/(2*N),  Fs1*(N-1)/(2*N),  Fs1*at)';

ok, i will rewrite it with an example:

If i sample 11 points (N=11) at the rate Fs=4 and fft it,
my first point on the frequency vector will be zero.
The next will be df=Fs/N=0.36364.
The 3rd will be 2*df=Fs/N=0.72727
...
...

In other words:

point k, frequency
point 1:   0  * 0.36364 =   0.00000
point 2:   1   *  0.36364 =   0.36364
point 2:   2   *  0.36364 =   0.72727
...
point 6:   5   *  0.36364 =   1.81818
point 7:   -1   *  0.36364 =   -0.36364
point 8:   -2   *  0.36364 =   -0.72727
...
point 11:   -5   *  0.36364 =   -1.81818

Points 1-6 are correct but 7-11 are backwards.  Point 7 should be -5 * 0.36364 and point 11 should be -1 * 0.36364.

Let's say you place the original points in order as follows:
[1 2 3 4 5 6 7 8 9 10 11]
When you do fftshift() on the fft output, the points get rearranged as follows:
[7 8 9 10 11 1 2 3 4 5 6]

So after this rearrangement, the frequencies must be in this order:
(-5:5)*0.36364

I hope this makes sense!

Again in other words, here is the code:

Fs=4                          % sampling rate
N=11                          % nr. of points
df=Fs/N                       % freq. increment per point
MaxFreq=df*(N-1)/2            % max measurable freq
freqIdx=[(0:N/2)'; -(1:N/2)'] % index vector to get the freqs
freq=freqIdx*df;              % frequency vector
[freqIdx freq]                % index beside frequency
fftshift(freq)                % shifted freq vector

All the best  ...

> Hi Brian,
>
>> freq=linspace(-Fs1*(N-1)/(2*N),  Fs1*(N-1)/(2*N),  Fs1*at)';
>
> ok, i will rewrite it with an example:
>
> If i sample 11 points (N=11) at the rate Fs=4 and fft it,
> my first point on the frequency vector will be zero.
> The next will be df=Fs/N=0.36364.
> The 3rd will be 2*df=Fs/N=0.72727
> ...
> ...
>
> In other words:
>
> point k, frequency
> point 1:   0  * 0.36364 =   0.00000
> point 2:   1   *  0.36364 =   0.36364
> point 2:   2   *  0.36364 =   0.72727
> ...
> point 6:   5   *  0.36364 =   1.81818
> point 7:   -1   *  0.36364 =   -0.36364
> point 8:   -2   *  0.36364 =   -0.72727
> ...
> point 11:   -5   *  0.36364 =   -1.81818
>
>
> Again in other words, here is the code:
>
>
> Fs=4                          % sampling rate
> N=11                          % nr. of points
> df=Fs/N                       % freq. increment per point
> MaxFreq=df*(N-1)/2            % max measurable freq
> freqIdx=[(0:N/2)'; -(1:N/2)'] % index vector to get the freqs
> freq=freqIdx*df;              % frequency vector
> [freqIdx freq]                % index beside frequency
> fftshift(freq)                % shifted freq vector
>
>
> All the best  ...
>
>
>> 2014-03-12 11:45 GMT+01:00 oxy <address@hidden>:
>>
>>> hi all,
>>>
>>> again in the question about FFT frequency axis.
>>>
>>> Pls confirm if I'm doing it right here:
>>>
>>> If I do fftshift(fft(time_domain_signal))
>>> of a certain time_domain_signal with even number of points,
>>> then the frequency axis becomes
>>>
>>> freq=linspace(-Fs1/2,  Fs1/2-Fs1/(Fs1*at),  Fs1*at)';
>>>
>>> where Fs= sampling rate, at=acquisition time.
>>>
>>> If however my time domain signal has an odd number of points, then
>>>
>>> freq=linspace(-Fs1/2,  Fs1/2,  Fs1*at)';
>>>
>>> Am id doing it right?
>>>
>>> We had a recent thread on that. See some msg below.
>>>
>>>
>>> oxy
>>>
>>>
>>> > % TTF Demo: this code shows how to calculate the
>>> > % frequency domain when doing FFT of td-signals.
>>> >
>>> > %--- parameter setup: same signal, 2 sampling rates ---
>>> > clear all
>>> > nu=10;            % signal frequency
>>> > at=5;             % acquisition time
>>> > T1=2              % signal decay constant in time domain
>>> > FsFactor1=4;      % sampling rate factor 1, Fs1/nu
>>> > FsFactor2=16;     % sampling rate factor 2, Fs2/nu
>>> >
>>> > %----------- calculating -------------------
>>> > Fs1 = FsFactor1*nu;            % Sampling rate 1
>>> > Fs2 = FsFactor2*nu;            % Sampling rate 2
>>> >
>>> > t1 = ((0:(at*Fs1-1))/Fs1)';    % Time vector 1
>>> > t2 = ((0:(at*Fs2-1))/Fs2)';    % Time vector 2
>>> >
>>> > tdsig1=exp(-i*2*pi*nu*t1).*exp(-t1./(T1)); % time domain signal 1
>>> > tdsig2=exp(-i*2*pi*nu*t2).*exp(-t2./(T1)); % time domain signal 2
>>> >
>>> > % I realized, the signal was appearing negative, thus minus...
>>> > freq1=linspace(-Fs1/2,Fs1/2-Fs1/(Fs1*at),Fs1*at)'; % freq. vector 1
>>> > freq2=linspace(-Fs2/2,Fs2/2-Fs2/(Fs2*at),Fs2*at)'; % freq. vector 2
>>> >
>>> > freqsig1=fftshift(fft(tdsig1));  % signal vector 1
>>> > freqsig2=fftshift(fft(tdsig2));  % signal vector 2
>>> >
>>> >
>>> > %----------- ploting -------------------
>>> > clf
>>> >   subplot(1,2,1)
>>> >   plot(freq1, freqsig1)
>>> >   axis([-20 0]), title('FsFactor=4')
>>> >
>>> >   subplot(1,2,2)
>>> >   plot(freq2, freqsig2)
>>> >   axis([-20  0]), title('FsFactor=16')
>>> > % --------- end ---------------------------
>>> >
>>> >>>================================================
>>> >>> Just one question more: you are calculating the frequency
>>> >>> vector as:
>>> >>>
>>> >>>          freq1=-linspace(-Fs1/2,   Fs1/2-1/at,   Fs1*at)';
>>> >>>
>>> >>> that makes it non-symetric. Isn't it better ...
>>> >>>
>>> >>>          freq1=-linspace(-Fs1/2+1/at,   Fs1/2-1/at,   Fs1*at)';
>>> >>>
>>> >>> ...? Then it's symetric. Hey, thanks a lot! Great help!!!!!
>>> >>
>>> >> That would be "better" if you really want to insist on having a
>>> >> symmetrical frequency vector at the expense of accuracy.
>>> >> Unfortunately, it's inaccurate.  For one thing, your symmetrical
>>> >> frequency vector does not contain a DC term if you have an even
>>> >> number
>>> >> of FFT points (which you should based on your examples).
>>> >>
>>> >> If you want the FFT value at the nth point to really correspond to
>>> >> freq1(n) you should use the slightly offset frequency vector that I
>>> >> provided.
>>> >>
>>> >> If it is more important to whatever you're doing that the frequency
>>> >> vector be exactly symmetrical and you can live with sacrificing the
>>> >> true 1-to-1 correspondence between frequency points and FFT points,
>>> >> you can use your symmetrical frequency vector.
>>> >>
>>> >>> oxy
>>> >>>
>>> >>> =================================================
>>> >>> On 2/19/14, oxy <address@hidden> wrote:
>>> >>>> On 2/12/14, Brian Kaczynski <address@hidden> wrote:
>>> >>>>> Hi Oxy,
>>> >>>>>
>>> >>>>> Since you're taking the FFT of a complex vector it will have
>>> >>>>> unique
>>> >>>>> terms
>>> >>>>> for negative vs. positive frequency.  You are correct that it
>>> >>>>> makes
>>> >>>>> more
>>> >>>>> sense to think of the FFT from -Fs/2 to +Fs/2.  It's just that the
>>> >>>>> Octave
>>> >>>>> fft function doesn't compute the bins in that order.  Due to
>>> aliasing,
>>> >>>>> any
>>> >>>>> frequency Fs/2 + df is equivalent to (-Fs/2 + df) so it's more a
>>> >>>>> matter
>>> >>>>> of
>>> >>>>> choice how you want to display the data.
>>> >>>>>
>>> >>>>> If you prefer plotting from -Fs/2 to +Fs/2 I would change these
>>> >>>>> two
>>> >>>>> lines
>>> >>>>> of your code as follows:
>>> >>>>>
>>> >>>>> freq=linspace(-Fs/2,Fs/2-Fs/(Fs*at),Fs*at)'
>>> >>>>> freqsig=fftshift(fft(tdsig);
>>> >>>>>
>>> >>>>> Let us know if that gives you what you want!
>>> >>>>>
>>> >>>>> -Brian
>>> >>>>>
>>> >>>>>
>>> >>>>> 2014-02-12 17:15 GMT+01:00 oxy <address@hidden>:
>>> >>>>>
>>> >>>>>> Hi Brian,
>>> >>>>>>
>>> >>>>>> > The frequency vector for FFT should go from DC as bin 1 to
>>> slightly
>>> >>>>>> > less
>>> >>>>>> > than Fs in the last bin (Fs - Fs/N where N is the number of FFT
>>> >>>>>> > points).
>>> >>>>>>
>>> >>>>>> If i understand u correctly, the frequency vector (freq) must be
>>> >>>>>> written
>>> >>>>>> as in this modified version of code below. However, according to
>>> >>>>>> Nyquist
>>> >>>>>> we
>>> >>>>>> cannot measure a frequency higher than Fs/2. Thus i do not see
>>> >>>>>> the meaning of plotting up to ~Fs.
>>> >>>>>>
>>> >>>>>> Also, if I rerun the code below changing FsFactor, I again see
>>> >>>>>> this
>>> >>>>>> dependency of the signal frequency on FsFactor.
>>> >>>>>> I cannot understand it. Looks like basics, yet not that obvious.
>>> >>>>>>
>>> >>>>>>   #---------- start code ------------
>>> >>>>>>   clear all
>>> >>>>>>   nu=10;                    % signal frequency
>>> >>>>>>   at=5;                       % acquisition time
>>> >>>>>>   T1=2                       % signal decay constant in time
>>> >>>>>> domain
>>> >>>>>>   FsFactor=16;             % the ratio (sampling rate)/(signal
>>> >>>>>> frequency), or Fs/nu
>>> >>>>>>
>>> >>>>>>   clf
>>> >>>>>>   Fs = FsFactor*nu;                             % Sampling rate
>>> >>>>>>   t = ((0:(at*Fs-1))/Fs)';                         % Time vector
>>> >>>>>>   % freq=linspace(-1,1,Fs*at)' * Fs/2;         % frequency
>>> >>>>>> vector,
>>> >>>>>> first
>>> >>>>>> version
>>> >>>>>>   freq=linspace(0,1,Fs*at)' * Fs- Fs/(Fs*at);  % frequency vector
>>> >>>>>> suggested by Brian
>>> >>>>>>   tdsig=exp(-i*2*pi*nu*t).*exp(-t./(T1));    % time domain signal
>>> >>>>>>   freqsig=fft(tdsig);                                % freq.
>>> >>>>>> domain
>>> >>>>>> signal
>>> >>>>>>   subplot(1,2,1)
>>> >>>>>>   plot(t,tdsig)
>>> >>>>>>   axis([ 0 0.4])       % zooming time domain to see that
>>> >>>>>> period=1/nu
>>> >>>>>>   subplot(1,2,2)
>>> >>>>>>   plot(freq, freqsig)
>>> >>>>>>   #---------- end code ------------
>>> >>>>>>
>>> >>>>>> thx guys ...
>>> >>>>>>
>>> >>>>>>
>>> >>>>>> > 2014-02-12 14:22 GMT+01:00 oxy <address@hidden>:
>>> >>>>>> >
>>> >>>>>> >> hey guys,
>>> >>>>>> >>
>>> >>>>>> >> the (simple) code bellow is how i ve learned to do FFT
>>> >>>>>> >> according
>>> >>>>>> >> to
>>> >>>>>> >> several docs online. However i observe a dependency of the
>>> >>>>>> >> signal
>>> >>>>>> >> frequency in the spectrum on the constant FsFactor. In other
>>> >>>>>> >> words,
>>> >>>>>> >> the signal frequency depends on the sampling rate. I should'nt
>>> be,
>>> >>>>>> >> right? So what is wrong here?
>>> >>>>>> >>
>>> >>>>>> >>   #---------- start code ------------
>>> >>>>>> >>   clear all
>>> >>>>>> >>   nu=10;                    % signal frequency
>>> >>>>>> >>   at=5;                       % acquisition time
>>> >>>>>> >>   T1=2                       % signal decay constant in time
>>> >>>>>> >> domain
>>> >>>>>> >>   FsFactor=8;             % the ratio (sampling rate)/(signal
>>> >>>>>> >> frequency), or Fs/nu
>>> >>>>>> >>
>>> >>>>>> >>   clf
>>> >>>>>> >>   Fs = FsFactor*nu;                             % Sampling
>>> >>>>>> >> rate
>>> >>>>>> >>   t = ((0:(at*Fs-1))/Fs)';                         % Time
>>> >>>>>> >> vector
>>> >>>>>> >>   freq=linspace(-1,1,Fs*at)' * Fs/2;         % frequency
>>> >>>>>> >> vector
>>> >>>>>> >>   tdsig=exp(-i*2*pi*nu*t).*exp(-t./(T1));    % time domain
>>> >>>>>> >> signal
>>> >>>>>> >>   freqsig=fft(tdsig);                                % freq.
>>> >>>>>> >> domain
>>> >>>>>> >> signal
>>> >>>>>> >>   subplot(1,2,1)
>>> >>>>>> >>   plot(t,tdsig)
>>> >>>>>> >>   axis([ 0 0.4])       % zooming time domain to see that
>>> >>>>>> >> period=1/nu
>>> >>>>>> >>   subplot(1,2,2)
>>> >>>>>> >>   plot(freq, freqsig)
>>> >>>>>> >>   #---------- end code ------------
>>> >>>>>> >>
>>> >>>>>> >> Do it yourself. Just rerun the code trying different values of
>>> >>>>>> >> FsFactor (eg: 2, 4, 8).
>>> >>>>>> >> Thx a lot for any hint!!!
>>> >>>>>> >>
>>> >>>>>> >> oxy
>>> >>>>>> >>
>>> >>>>>> >> ps: cross post
>>> >>>>>> >>
>>> >>>>>>
>>> >>>>>> >> _______________________________________________
>>> >>>>>> >> Help-octave mailing list
>>> >>>>>> >> https://mailman.cae.wisc.edu/listinfo/help-octave
>>> >>>>>> >>
>>> >>>>>> >
>>> >>>>>> _______________________________________________
>>> >>>>>> Help-octave mailing list
>>> >>>>>> https://mailman.cae.wisc.edu/listinfo/help-octave
>>> >>>>>>
>>> >>>>>
>>> >>>>
>>> >>>
>>> >>
>>> >
>>>
>>
>
_______________________________________________
Help-octave mailing list