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

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?

[Juan Manriquez]

Hello I did this experiment while to beat with this, I hope this to be usuful.

Implementation using example of http://www.robots.ox.ac.uk/~sjrob/Teaching/SP/l7.pdf

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.

Answer 1

[Teja Muppirala]

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'})

[Jeff Fischer]

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.

[Teja Muppirala]

Read the appended portion above.

[Indiase Straal]

Very insightful answer, thanks!

Answer 2

[Andrew Newell]

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.

[Jeff Fischer]

Try doing it with a longer FFT y=fft(x,2000)

[Andrew Newell]

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.

Answer 3

[Daniel Shub]

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

[Jeff Fischer]

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

Answer 4

[Walter Roberson]

Round off error.