MATLAB Answers

0

FFT result does not jive with theory for basic sine and cosine

Asked by Jeff Fischer on 15 May 2011
Latest activity Commented on by Kiran Anginthaya on 27 Nov 2018
If I create a simple sin(2*pi*.1*t) or same for a cos, and let t=[-50:50], and transform it, I get spikes in both the real and imaginary outputs. But Fourier transform theory says I should get only real components for Cosine and only Imaginary components for sine. Why is Matlab different?

  1 Comment

Hello I did this experiment while to beat with this, I hope this to be usuful.
clear all;
%
N=8;
n=zeros(N,1);
armonicos=zeros(size(n));
espectro=zeros(size(n));
%
for index=1:N
armonicos(index)=index-1;
end
%Durante el n se crea un array con los valores de la intensidad
%de la se?al evaluada en los puntos definidos al particionar el periodo
%de n.
for index=1:N
sample=index-1;
n(index)=signal_sampling(sample);
end
%Se transforma cada valor de intensidad muestreada de la se?al original
%en un multiplo del armonico fundamental (si este se encuentra en la)
%se?al original y se crea la tabla con los valores del espectro (modulo
%del vector complejo) obtenido para cada muestra.
for index=1:N
espectro(index)=abs(transformar_muestra(armonicos(index),n));
end
%
%------Graficacion del espectro en el dominio de frecuencias---------------
figure(1);stem(armonicos,espectro,'fill','--');
%plot(armonicos,espectro)
%--------------------------------------------------------------------------
%
%
%
%
%
%
%
%
%
%Aqui se define la funcion de la se?al de entrada
function sampling=signal_sampling(sample)
%sampling=1+cos((pi/2)*sample)+(1/2)*sin((6*pi/8)*sample);
sampling=cos(sample*pi/2);
%sampling=cos(sample);
%sampling=sample;
%sampling=5+2*cos(2*pi*sample-pi/2)+3*cos(4*pi*sample);
end
%Aqui se realiza la transformacion de cada muestra al dominio de
%frecuencias.
function resultado=transformar_muestra(armonico,n)
N=size(n,1);
k=1;
tetha=2*pi/N;
fasor=exp(-i*(armonico)*k*tetha);
e=zeros(size(n));
for index=1:N
multiplo=index-1;
e(index)=fasor^multiplo;
end
resultado=round(n.'*e)/N;
end
In fact I had the same problem but finally I achieve to get the same with matlab, but really I don't understand qhy using the double of the period is that this works
next is the code for doing this one using fft
x=(0:1:8);
x1=x(1,1:8);
y=1*cos((pi/2)*x1);
f=fft(y);
k=(0:1:7);
figure(3),stem(k,abs(f)/8);
it is everything I have.

Sign in to comment.

4 Answers

Answer by Teja Muppirala
on 16 May 2011
 Accepted Answer

There is absolutely nothing strange here.The discrete fourier transform assumes PERIODIC SIGNALS.
In your case:
t=[-50:50];
y = sin(2*pi*.1*t) ;
plot(t,y);
Y = fft(y);
you are not really taking the fourier transform of a single sine wave. The first element y(1) is zero, and the last element y(101) is zero, and therefore the periodic extension [ y y y y ...] has two consecutive zeros in it. When you take the fourier transform of that signal, that causes a whole bunch of other frequency components to come in.
If you do this however,
t=[-50:49];
y = sin(2*pi*.1*t) ;
plot(t,y);
Y = fft(y);
Then [y y y y y ...] is indeed a pure single periodic sinusoid, and you see that you get exactly the answer you were hoping for (a single frequency component with sine as purely imaginary, and cosine as purely real). Another example:
figure,
t=[0:0.01:1];
y = cos(2*pi*10*t) ;
Y = fft(y);
plot(real(Y));
hold all
t=[0:0.01:0.99];
y = cos(2*pi*10*t) ;
Y = fft(y);
plot(real(Y),'r');
legend({'An "almost" periodic cosine wave' 'A periodic cosine wave'})
EDIT:
Again, the answer lies in the fact that the discrete fourier transform assumes periodic signals.
Question: When you type this,
t=[-50:49];
y = sin(2*pi*.1*t) ;
Y=fft(y,2000);
What are you really taking the discrete fourier transform of? Answer: You are taking the fourier transform of this infinitely extended periodic signal [y zeros(1,1900) y zeros(1,1900) y zeros(1,1900) ...] .
Plot it:
plot([y zeros(1,1900) y zeros(1,1900) y zeros(1,1900)]);
Is this a sine wave? Clearly it is not. This curve has a bunch of components in it, some of which are real and some of which are imaginary. And when you plot the FFT, you are infact seeing all of those components.
Compare again:
figure
t=[0:0.01:19.99];
y = sin(2*pi*10*t);
Y = fft(y,2000);
plot(real(Y));
hold all
t=[0:0.01:0.99];
y = sin(2*pi*10*t) ;
Y = fft(y,2000);
plot(real(Y));
hold all
t=[0:0.01:19.99];
y = sin(2*pi*10*t);
y(101:end) = 0;
Y = fft(y);
plot(real(Y),'r:');
legend({'Real part of FFT of a sine wave' 'Real part of FFT of something that''s not a sine wave', 'The same Real part FFT of something that''s not a sine wave'})

  3 Comments

Thanks for trying but that isn't quite the answer.
Do the same thing with Y=fft(y,2000) and you will see just as significant a impulse in the real as the imaginary plane.

Sign in to comment.


Answer by Andrew Newell
on 15 May 2011

I don't think it is a round off error. If you run this code:
n = 1000;
t=linspace(0,2,n);
x=sin(2*pi*t);
y = fft(x);
x2 = ifft(y);
you'll find that x and x2 agree within machine precision. I think the spikes occur because this is a discrete Fourier transform, not the integral.
EDIT: To go a bit further: the Fourier integral transforms are F(cos(2*pi*k0*t)) = 0.5*(Dirac(k+k0)+Dirac(k-k_0)) F(sin(2*pi*k0*t)) = 0.5*1i*(Dirac(k+k0)-Dirac(k-k_0)). Now if we choose some random value for k0, and do the discrete Fourier transform of the cosine over [-k0,k0],
n = 1000;
k0 = rand(1);
t=linspace(-1,1,n);
x=cos(2*pi*k0*t);
y = fft(x);
plot(t/k0,real(y),t/k0,imag(y))
We find two positive spikes in the real component at t=-k0 and t=k0 and a much smaller one for the imaginary component. If we do the same for the sin,
t=linspace(-1,1,n);
x=sin(2*pi*k0*t);
y = fft(x);
figure
plot(t/k0,real(y),t/k0,imag(y))
we get a negative and positive spike in the imaginary component and much smaller spikes in the real component - much as predicted for the integrals.

  2 Comments

Try doing it with a longer FFT y=fft(x,2000)
That removes the peak on the right. I think, if anything, this supports my main point, that discrete Fourier transforms are similar to continuous Fourier, but also have some systematic differences.

Sign in to comment.


Answer by Daniel Shub
on 15 May 2011

Steve Eddins recently blogged all about FFT in MATLAB. http://blogs.mathworks.com/steve/category/fourier-transforms/

  1 Comment

Read through the blog, thanks. But he didn't seem to cover this.

Sign in to comment.


Answer by Walter Roberson
on 15 May 2011

Round off error.

  0 Comments

Sign in to comment.