FFT result does not jive with theory for basic sine and cosine
18 views (last 30 days)
Show older comments
Jeff Fischer
on 15 May 2011
Commented: Indiase Straal
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
Juan Manriquez
on 1 Jun 2017
Edited: Walter Roberson
on 1 Jun 2017
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.
Accepted Answer
Teja Muppirala
on 16 May 2011
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
More Answers (3)
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
Andrew Newell
on 16 May 2011
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.
Daniel Shub
on 15 May 2011
Steve Eddins recently blogged all about FFT in MATLAB. http://blogs.mathworks.com/steve/category/fourier-transforms/
See Also
Categories
Find more on Spectral Measurements in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!