MATLAB fft vs Mathematical fourier

Question

I am using the following code to generate a fft and mathematical Fourier transform of a signal. I want to then mathematically recreate the original signal of the fft. This works on the mathematical signal but not on the fft since it is a Discrete Transform. Does anyone know what change I can make to my inverse transform equation that will make it work for fft?

clear all; clc;
N = 1024;
N2 = 1023;
SNR = -10;
fs = 1024;
Ts = 1/fs;
t = (0:(N-1))*Ts;
x = 0.5*sawtooth(2*2*pi*t);
x1 = fft(x); 
Magnitude1 = abs(x1);
Phase1 = angle(x1)*360/(2*pi);

for m = 1:1024
   f(m) = m;                % Sinusoidal frequencies
   a = (2/N)*sum(x.*cos(2*pi*f(m)*t));      % Cosine coeff. 
   b = (2/N)*sum(x.*sin(2*pi*f(m)*t));       % Sine coeff
   Magnitude(m) = sqrt(a^2 + b^2);                % Magnitude spectrum
   Phase(m) = -atan2(b,a);                   % Phase spectrum
end

subplot(2,1,1);
plot(f,Magnitude1./512);   % Plot magnitude spectrum
       ......Labels and title.......
subplot(2,1,2);
plot(f,Magnitude,'k');    % Plot phase spectrum
ylabel('Phase (deg)','FontSize',14);
pause();

x2 = zeros(1,1024);         % Waveform vector 
for m = 1:24
    f(m) = m;               % Sinusoidal frequencies
    x2 = (1/m)*(x2 + Magnitude1(m)*cos(2*pi*f(m)*t + Phase1(m)));  
end
x3 = zeros(1,1024);         % Waveform vector
for m = 1:24
    f(m) = m;               % Sinusoidal frequencies
    x3 = (x3 + Magnitude(m)*cos(2*pi*f(m)*t + Phase(m)));  
end
plot(t,x,'--k'); hold on;
plot(t,x2,'k');
plot(t,x3,'b');```

Thales Thales · Accepted Answer · 2019-10-07T16:57:36

There are a few comments about the Fourier Transform, and I hope I can explain everything for you. Also, I don't know what you mean by "Mathematical Fourier transform", as none of the expressions in your code is resembles the Fourier series of the sawtooth wave.

To understand exactly what the fft function does, we can do things step by step. First, following your code, we create and plot one period of the sawtooth wave.

n = 1024;
fs = 1024;
dt = 1/fs;
t = (0:(n-1))*dt;
x = 0.5*sawtooth(2*pi*t);

figure; plot(t,x); xlabel('t [s]'); ylabel('x');

We can now calculate a few things. First, the Nyquist frequency, the maximum detectable frequency from the samples.

f_max = 0.5*fs

 f_max =

    512

Also, the minimum detectable frequency,

f_min = 1/t(end)

 f_min =

    1.000977517106549

Calculate now the discrete Fourier transform with MATLAB function:

X = fft(x)/n;

This function obtains the complex coefficients of each term of the discrete Fourier transform. Notice it calculates the coefficients using the exp notation, not in terms of sines and cosines. The division by n is to guarantee that the first coefficient is equal to the arithmetic mean of the samples

If you want to plot the magnitude/phase of the transformed signal, you can type:

f = linspace(f_min,f_max,n/2); % frequency vector
a0 = X(1); % constant amplitude
X(1)=[]; % we don't have to plot the first component, as it is the constant amplitude term
XP = X(1:n/2); % we get only the first half of the array, as the second half is the reflection along the y-axis

figure
subplot(2,1,1)
plot(f,abs(XP)); ylabel('Amplitude');
subplot(2,1,2)
plot(f,angle(XP)); ylabel('Phase');
xlabel('Frequency [Hz]')

What does this plot means? It shows in a figure the amplitude and phase of the complex coefficients of the terms in the Fourier series that represent the original signal (the sawtooth wave). You can use this coefficients to obtain the signal approximation in terms of a (truncated) Fourier series. Of course, to do that, we need the whole transform (not only the first half, as it is usual to plot it).

X = fft(x)/n;
amplitude = abs(X);
phase = angle(X);
f = fs*[(0:(n/2)-1)/n (-n/2:-1)/n]; % frequency vector with all components

% we calculate the value of x for each time step
for j=1:n
    x_approx(j) = 0;
    for k=1:n % summation done using a for
        x_approx(j) = x_approx(j)+X(k)*exp(2*pi*1i/n*(j-1)*(k-1));
    end
    x_approx(j) = x_approx(j);
end

Notice: The code above is for clarification and does not intend to be well coded. The summation can be done in MATLAB in a much better way than using a for loop, and some warnings will pop up in the code, warning the user to preallocate each variable for speed.

The above code calculates the x(ti) for each time ti, using the terms of the truncated Fourier series. If we plot both the original signal and the approximated one, we get:

figure
plot(t,x,t,x_approx)
legend('original signal','signal from fft','location','best')

The original signal and the approximated one are nearly equal. As a matter of fact,

norm(x-x_approx)

 ans =

      1.997566360514140e-12

Is almost zero, but not exactly zero. Also, the plot above will issue a warning, due to the use of complex coefficients when calculating the approximated signal:

Warning: Imaginary parts of complex X and/or Y arguments ignored

But you can check that the imaginary term is very close to zero. It is not exactly zero due to roundoff errors in the computations.

norm(imag(x_approx))

 ans =

      1.402648396024229e-12

Notice in the codes above how to interpret and use the results from the fft function and how they are represented in the exp form, not on terms of sines and cosines, as you coded.

MATLAB fft vs Mathematical fourier

1 Answers