DFT in Matlab without built-in function

The discrete Fourier transform (DFT) is a powerful tool for analyzing the frequency content of digital signals. It allows us to transform a sequence of N complex numbers into a sequence of N complex numbers that represent the signal's frequency components. 

Matlab has built-in function called fft() to calculate DFT. Learning to code without using inbuilt functions can be an important skill for several reasons. One of the reason is understanding the underlying principles: When we use inbuilt or built-in functions, we may not always be aware of the underlying algorithms and mathematical concepts that are being used to perform the calculations. However, when we code from scratch, we have to understand the principles involved in the calculations, which can deepen our understanding of the concepts and lead to better problem-solving skills.

In this digital signal processing with Matlab tutorial, we will explore how to implement the DFT in MATLAB without using inbuilt functions, that is, calculate DFT in matlab without using fft function.

Overview of DFT

The DFT of a sequence x[n] of length N is defined as follows:

\( X[k] = \sum_{n=0}^{ N-1}  x[n] e^{(\frac{-j  2 \pi }{N}nk)} \) ---->(1)

for \(0 \leq k \leq N-1\) 

where X[k] represents the complex amplitude of the k-th frequency component, and j is the imaginary unit.

The above equation can be implemented in MATLAB as follows:

N = length(x); % Length of input sequence
X = zeros(1,N); % Initialize output sequence

for k = 0:N-1 % Loop over all frequency components
    for n = 0:N-1 % Loop over all time-domain samples
        X(k+1) = X(k+1) + x(n+1)*exp(-1i*2*pi*n*k/N); % Calculate DFT
    end
end

In this code, N is the length of the input sequence x, and X is an array of zeros that will be used to store the DFT of x. We loop over all frequency components k and time-domain samples n, and calculate the DFT using the formula above.

Magnitude and Phase

Once we have calculated the DFT of a signal, we can obtain its magnitude and phase by taking the absolute value and angle of each frequency component, respectively.

The magnitude of the k-th frequency component is given by:

mag_X(k) = sqrt(real(X[k])^2 + imag(X[k])^2)

and the phase is given by:

DFT Matlab Code without Built-In Function

Below is an example MATLAB code to compute the Discrete Fourier Transform (DFT) of a given signal x without using any built-in functions and using the formula mag_X = sqrt(real(X[k])^2 + imag(X[k])^2) to compute the magnitude and phi_X = atan2(imag(X[k]), real(X[k])) to compute the phase:

% Input signal
fs = 100; % Sampling frequency
t = 0:1/fs:1-1/fs; % Time vector
x = sin(2*pi*15*t) + sin(2*pi*40*t); % Signal

% Compute DFT
N = length(x); % Length of signal
X = zeros(1,N); % Initialize DFT coefficients
for k = 1:N
    for n = 1:N
        X(k) = X(k) + x(n)*exp(-1i*2*pi*(k-1)*(n-1)/N);
    end
end

% Compute magnitude and phase
mag_X = zeros(1,N); % Initialize magnitude
phi_X = zeros(1,N); % Initialize phase
for k = 1:N
    mag_X(k) = sqrt(real(X(k))^2 + imag(X(k))^2); % Compute magnitude
    phi_X(k) = atan2(imag(X(k)), real(X(k))); % Compute phase
end

% Plot results
f = (0:N-1)*fs/N; % Frequency vector
subplot(2,1,1);
plot(f, mag_X);
xlabel('Frequency (Hz)');
ylabel('Magnitude');
title('Magnitude Spectrum');
subplot(2,1,2);
plot(f, phi_X);
xlabel('Frequency (Hz)');
ylabel('Phase (rad)');
title('Phase Spectrum');

In this code, we first define the input signal x as a sum of two sinusoidal components with frequencies of 15 Hz and 40 Hz, respectively. We then compute the DFT coefficients X using two nested for loops that implement the formula for the DFT.

Next, we compute the magnitude and phase of each frequency component in the DFT using the formulas |X[k]| = sqrt(real(X[k])^2 + imag(X[k])^2) and phi[k] = atan2(imag(X[k]), real(X[k])), respectively. We store the results in the arrays mag_X and phi_X.

Finally, we plot the magnitude and phase spectra using the plot() function in MATLAB. The frequency axis is defined using the sampling frequency fs and the length of the signal N.

The following is DFT plot, that is the magnitude response and phase response.

In this way we can perform DFT calculation in Matlab without using in built functions such as fft(), abs(), unwrap() or angle() functions. The reason for not using built-in function is to understand the underlying mathematics. Other reasons and advantages are:

• Flexibility and customization: Inbuilt functions may not always have the specific functionality that we need, or we may need to modify them to fit our specific requirements. When we know how to code from scratch, we have the flexibility to customize the code to our specific needs and make modifications as necessary.
• Debugging and troubleshooting: When we use inbuilt functions, it can be difficult to debug and troubleshoot errors that may arise. However, when we code from scratch, we have more control over the code and can easily identify and fix errors that may arise.
• Portability: Inbuilt functions may not always be available or compatible with different versions of software or programming languages. When we know how to code from scratch, we can write code that is portable and can be used across different platforms and programming languages.

Overall, learning to code without using inbuilt functions can improve our problem-solving skills, deepen our understanding of the underlying concepts, and give us greater flexibility and control over our code.

References:

[1] What is Discrete Fourier Transform(DFT)

[2] DFT in Matlab without built-in function