Thread Subject: FFT acting strangely...

I am trying to perform an FFT on some data. My data is:

times = [1:1000];
accelArray = 0.7*sin(50*times);

Then I do:

fftAccel = (fft(accelArray))/length(accelArray);
fftFreqs = Fs/2*linspace(0,1,length(accelArray)/2);

Then, when plotting, I do this:

plot(fftFreqs,2*abs(fftAccel(1:length(accelArray)/2)));

Now, I'm expecting to get a peak at 50Hz, but instead I get
a peak at 42.8Hz. I've tried to put sin(2*pi*50*times), but
that just gives me numbers that are close to 10^-11.

Any ideas?

Mike

On 29 Feb, 20:56, "Michael Stachowsky" <mstachow...@gmail.com> wrote:
> I am trying to perform an FFT on some data. =A0My data is:
>
> times =3D [1:1000];
> accelArray =3D 0.7*sin(50*times);
>
> Then I do:
>
> fftAccel =3D (fft(accelArray))/length(accelArray);
> fftFreqs =3D Fs/2*linspace(0,1,length(accelArray)/2);
>
> Then, when plotting, I do this:
>
> plot(fftFreqs,2*abs(fftAccel(1:length(accelArray)/2)));
>
> Now, I'm expecting to get a peak at 50Hz, but instead I get
> a peak at 42.8Hz. =A0I've tried to put sin(2*pi*50*times), but
> that just gives me numbers that are close to 10^-11.

I tried to run your code, but got an error message that Fs is not
specified. You need to specify it and make sure that the vector
times is set up as

Fs =3D <????>
times =3D (0:1000)/Fs;

Rune

Ah, sorry. Fs = 1000;

I discovered my problem: it was one of aliasing. I had
mistyped the example data.

But I came upon another issue. This may stem from my lack
of familiarity with FFT, but, how do I make the proper array
of frequencies for my analysis?

This is going to be used in an experiment with reasonably
strict accuracy controls (at least +/- 0.1Hz), and my FFT
analysis gave me 120.12Hz when I gave it a 120Hz sine wave...

Mike

Rune Allnor <allnor@tele.ntnu.no> wrote in message
<f187d8bb-5e73-491c-94ae-56d64264cfbe@p25g2000hsf.googlegroups.com>...
> On 29 Feb, 20:56, "Michael Stachowsky"
<mstachow...@gmail.com> wrote:
> > I am trying to perform an FFT on some data. =A0My data is:
> >
> > times =3D [1:1000];
> > accelArray =3D 0.7*sin(50*times);
> >
> > Then I do:
> >
> > fftAccel =3D (fft(accelArray))/length(accelArray);
> > fftFreqs =3D Fs/2*linspace(0,1,length(accelArray)/2);
> >
> > Then, when plotting, I do this:
> >
> > plot(fftFreqs,2*abs(fftAccel(1:length(accelArray)/2)));
> >
> > Now, I'm expecting to get a peak at 50Hz, but instead I get
> > a peak at 42.8Hz. =A0I've tried to put
sin(2*pi*50*times), but
> > that just gives me numbers that are close to 10^-11.
>
> I tried to run your code, but got an error message that Fs
is not
> specified. You need to specify it and make sure that the
vector
> times is set up as
>
> Fs =3D <????>
> times =3D (0:1000)/Fs;
>
> Rune

Top-posting repaired.

On Feb 29, 3:51 pm, "Michael Stachowsky" <mstachow...@gmail.com>
wrote:
> Rune Allnor <all...@tele.ntnu.no> wrote in message
>
> <f187d8bb-5e73-491c-94ae-56d64264c...@p25g2000hsf.googlegroups.com>...
>
> > On 29 Feb, 20:56, "Michael Stachowsky"
> <mstachow...@gmail.com> wrote:
> > > I am trying to perform an FFT on some data. =A0My data is:
>
> > > times =3D [1:1000];
> > > accelArray =3D 0.7*sin(50*times);
>
> > > Then I do:
>
> > > fftAccel =3D (fft(accelArray))/length(accelArray);
> > > fftFreqs =3D Fs/2*linspace(0,1,length(accelArray)/2);
>
> > > Then, when plotting, I do this:
>
> > > plot(fftFreqs,2*abs(fftAccel(1:length(accelArray)/2)));
>
> > > Now, I'm expecting to get a peak at 50Hz, but instead I get
> > > a peak at 42.8Hz. =A0I've tried to put
>
> sin(2*pi*50*times), but
>
> > > that just gives me numbers that are close to 10^-11.

The reason this happened is that the way you
defined "times" was as a set of integers. So
you are calculating sin(2*pi*50*integer), which
is the sine of an integer multiple of 2*pi, which is
0.

> > I tried to run your code, but got an error message that Fs
> is not
> > specified. You need to specify it and make sure that the
> vector
> > times is set up as
>
> > Fs = <????>
> > times = (0:1000)/Fs;
>
> Ah, sorry. Fs = 1000;

sin(2*pi*f*times) is correct for a sine wave of frequency
f, but if you want to sample at rate Fs, then your times
should be multiples of 1/Fs as Rune suggested.

If you want 1000 samples, then you should use
  times = (0:999)/Fs;

In general, if you want N samples at rate Fs, your times
are
   times = (0:N-1)/Fs;

> I discovered my problem: it was one of aliasing. I had
> mistyped the example data.
>
> But I came upon another issue. This may stem from my lack
> of familiarity with FFT, but, how do I make the proper array
> of frequencies for my analysis?

The frequency spacing is 1/(time of sample). If you have
N samples at rate Fs, then it is 1/(N/Fs) = Fs/N

Thus, your frequency array should be (0:N-1) * (Fs/N)

> This is going to be used in an experiment with reasonably
> strict accuracy controls (at least +/- 0.1Hz), and my FFT
> analysis gave me 120.12Hz when I gave it a 120Hz sine wave...

My guess is that this is because you used a 1001-point
time sample instead of a 1000-point time sample, since
your frequency seems to be off by 1 part in 1000.

                 - Randy

On 29 Feb, 21:51, "Michael Stachowsky" <mstachow...@gmail.com> wrote:
> Ah, sorry. =A0Fs =3D 1000;
>
> I discovered my problem: it was one of aliasing. =A0I had
> mistyped the example data.
>
> But I came upon another issue. =A0This may stem from my lack
> of familiarity with FFT, but, how do I make the proper array
> of frequencies for my analysis?

You don't. The array of frequencies is given implicitly as

fv =3D (0:N-1)/N*Fs

where N is the length of the data sequence fed to the FFT
and Fs is the sampling frequency.

> This is going to be used in an experiment with reasonably
> strict accuracy controls (at least +/- 0.1Hz), and my FFT
> analysis gave me 120.12Hz when I gave it a 120Hz sine wave...

You have to determine a largest acceptable bin width df
and compute the minimum sequence length as

N >=3D Fs/df

For Fs =3D 1000 Hz and df =3D 0.1 Hz one gets

N >=3D 1000/0.1 =3D 10000

mening that you need at least 10000 samples in the
signal vector to guarantee that you end up within
0.1 Hz of the desired signal.

Rune

On Feb 29, 4:16 pm, Rune Allnor <all...@tele.ntnu.no> wrote:
> You don't. The array of frequencies is given implicitly as
>
> fv = (0:N-1)/N*Fs
>
> where N is the length of the data sequence fed to the FFT
> and Fs is the sampling frequency.

Both you and Randy gave this answer, which is the "obvious" choice if
you look at the DFT definition, but in most circumstances this is a
very bad definition for the array of frequencies.

Because of aliasing, the frequencies are ambiguous. However, if you
assume that your data are band-limited up to the Nyquist frequency,
then the correct array of frequencies (assuming even N) is:

fv = [0:N/2,-N/2+1:-1]/N*Fs

Another way of looking at it is that, by choosing the frequency
aliasing (adding any multiple of N to any of the integers in the [...]
array), you are specifying a trigonometric interpolation in between
your sampling points. The frequencies I gave above correspond to the
minimal-slope interpolation (except for the N/2 Nyquist frequency,
which has to be handled specially if you really care about it), which
is usually the best choice.

(If the data are real, then the second half of the spectrum is
redundant and you might as well only plot for [0:N/2].)

Regards,
Steven G. Johnson

On Feb 29, 3:51 pm, "Michael Stachowsky" <mstachow...@gmail.com>
wrote:
> Ah, sorry. Fs = 1000;
>
> I discovered my problem: it was one of aliasing. I had
> mistyped the example data.
>
> But I came upon another issue. This may stem from my lack
> of familiarity with FFT, but, how do I make the proper array
> of frequencies for my analysis?

It depends on whether you have an odd or an even number of points.

dt = 1/Fs % Time sample spacing
T = N*dt % FFT imposed periodicity
t = 0: dt: T-dt; % t = (0:N-1)*dt;

df = 1/T % Frequency sample spacing
Fs = N*df % FFT imposed periodicity
f = 0: df: Fs-df % f = (0:N-1)*df;

f0 = 50;
x = 0.7*sin(2*pi*50*t);
X = fft(x); % Defined on unipolar frequency interval [0,Fs-df]
Xb = fftshift(X); % Defined on bipolar frequency interval [-?,+?]

The bipolar frequency range depends on whether N is odd or even

Q = ceil((N+1)/2)
P = N+1-Q % P = floor((N+1)/2)
fQ = (Q-1)*df % Maximum negative frequency
fP = (P-1)*df % Maximum positive frequency

Then

f = 0: df: fP+fQ; % Unipolar, used with X
fb = f - fQ;
fb = -fQ: df : fP; % Bipolar, used with Xb

Notice that when N is even

fNyquist = fQ > fP

So MATLAB's fftshift puts the Nyquist frequency in the
negative part of the spectrum.


> This is going to be used in an experiment with reasonably
> strict accuracy controls (at least +/- 0.1Hz), and my FFT
> analysis gave me 120.12Hz when I gave it a 120Hz sine wave...

Since

df = 1/T = 1/(N*dt) = Fs/N,

There are 2 ways to decrease df.


Hope this helps.

Greg

On Feb 29, 4:54=A0pm, "Steven G. Johnson" <stev...@alum.mit.edu> wrote:
> On Feb 29, 4:16 pm, Rune Allnor <all...@tele.ntnu.no> wrote:
>
> > You don't. The array of frequencies is given implicitly as
>
> > fv =3D (0:N-1)/N*Fs
>
> > where N is the length of the data sequence fed to the FFT
> > and Fs is the sampling frequency.
>
> Both you and Randy gave this answer, which is the "obvious" choice if
> you look at the DFT definition, but in most circumstances this is a
> very bad definition for the array of frequencies.
>
> Because of aliasing, the frequencies are ambiguous. =A0However, if you
> assume that your data are band-limited up to the Nyquist frequency,
> then the correct array of frequencies (assuming even N) is:
>
> fv =3D [0:N/2,-N/2+1:-1]/N*Fs

It can be very confusing to try to understand this without

1. The mention of fftshift
2. The explicit restriction to N even (The OP had N =3D 1001)

> Another way of looking at it is that, by choosing the frequency
> aliasing (adding any multiple of N to any of the integers in the [...]
> array), you are specifying a trigonometric interpolation in between
> your sampling points. =A0
>The frequencies I gave above correspond to the
> minimal-slope interpolation (except for the N/2 Nyquist frequency,
> which has to be handled specially if you really care about it), which
> is usually the best choice.

Isn't the interpolation in terms of sinc functions?

> (If the data are real, then the second half of the spectrum is
> redundant and you might as well only plot for [0:N/2].)

Clarification:

Q =3D ceil(N+1)/2), fQ =3D (Q-1)*df

"First half" of spectrum is

X(1:Q) , f(1:Q) =3D [0:Q-1]*df =3D 0: df: fQ

where

Q =3D (N+1)/2 for N odd
Q =3D 1+N/2 for N even

Hope this helps.

Greg

On Mar 1, 7:42 pm, Greg Heath <he...@alumni.brown.edu> wrote:
> Isn't the interpolation in terms of sinc functions?

One way to look at the interpolation is as a convolution of the
original data with sinc functions. However, it is completely
equivalent to look at the interpolation as simple sines and cosines
(complex exponentials) multiplied by the Fourier amplitudes -- it
depends on whether you want to express the interpolation (in time
domain) in terms of the time- or frequency-domain coefficients. The
expression in terms of the frequency-domain coefficients is much
easier to analyze (being a simple sum and not a convolution, and being
in terms of orthogonal functions), and corresponds to a trigonometric
interpolation (a fit of the original data to sines and cosine, which
passes through all the original data points).

(Trigonometric interpolation, in fact, was the original context in
which both the DFT and the most common FFT algorithm were proposed by
Gauss---he was fitting asteroid observations to sines and cosines in
order to interpolate, since his data were periodic.)

Thanks for clarifying my answer regarding zero-based vs. one-based
indices; and I hadn't noticed that the OP's N was odd.

Steven

On Mar 3, 11:45 am, "Steven G. Johnson" <stev...@alum.mit.edu> wrote:
> expression in terms of the frequency-domain coefficients is much
> easier to analyze (being a simple sum and not a convolution,

Well, I should correct myself that the DFT can be written as a
convolution, and the direct expression doesn't have any fewer
operations than a convolution, so this isn't realy an advantage of one
expression over the other. But the fact that the trigonometric-
interpolation viewpoint is in terms of orthogonal functions certainly
makes it easier to analyzer for some purposes, as is the fact that
they are eigenfunctions of differentiation; both of these make it easy
to show which choice of aliasing leads to minimal rms slope, which
seems harder to do in terms of the sinc-convolution picture.