discrete fourier transform code in matlab

[Saad Muftah]

Hi everyone

I am trying to do matlab code to calculate the real and imaginary values of fundamental frequency sine wave by discrete fourier transform for one period of 16 sample. And i reached a dead end, the code that i used or in other word the when i was implementing the equation of fourier transform i didn't get the right answer. for example if i enter 5sin(377*t), and calculate the real and imaginary value give a strange numbers such as 26.

Please help me guys, if some one can help me.

Many thanks,


Smuftahndi@gmail.com

 

[peanutnjelly]

Here you have DFT. This code is for learning purpose or small signal size, because it gets really slow with huge signal size N (because its DFT).
And for DFT signal size is not necessary 2^N, like for FFT.

Code is for scilab(similar).

function y = my_DFT(sig,N)
    for k=1 : N-1
        re_sum =0
        im_sum =0
        for n=1:N-1
            cosf = cos((2*%pi*n*k)/N)
            sinf = sin((2*%pi*n*k)/N)
            re_sum = re_sum + real(sig(n))*cosf-imag(sig(n))*sinf
            im_sum = im_sum + real(sig(n))*sinf+imag(sig(n))*cosf
        end
        y(k)= re_sum + (%i * im_sum) 
    end
endfunction

scilab code to test it:

t=0:0.001:0.1
signal = sin(2*%pi*100*t)
plot(abs(my_DFT(signal,100)))

 

[Saad Muftah]

Thanks for your reply, but it didn't work for me, i don't know if you could upload it.

 

[weetabixharry]

Hi everyone

I am trying to do matlab code to calculate the real and imaginary values of fundamental frequency sine wave by discrete fourier transform for one period of 16 sample.

It's easy to make mistakes with this stuff if you're a beginner. I have written this code that should help you to understand what you need to do:

close all; clear all; clc;

% ============= (1) Time-domain Signal ============= %

% Frequency of the sine in Hertz
f_sine = 377;

% Period of the sine
T_sine = 1/f_sine;

% Create "Nsamps" equally-spaced points in time from 0 to T_sine
Nsamps = 16;
t = (0:Nsamps-1)*T_sine/Nsamps;

% Period between each sample of t
T_samp = t(2)-t(1);
f_samp = 1/T_samp; % Sampling frequency

% Create the sine signal at time points t
x = sin(2*pi*f_sine*t);

% Plot time-domain signal
figure();
plot(t,x);
title('Time-domain signal, x(t)');
xlabel('Time');
ylabel('Amplitude');
grid on; axis tight;

% ========== (2) Frequency-domain Signal =========== %

% Frequency difference between each frequency bin in spectrum
df = f_samp/Nsamps;

% All frequencies to plot in spectrum
f = df*(0:Nsamps-1);

% Discrete Fourier Transform (using FFT method)
Y = fft(x);

% Plot frequency spectrum...
figure();
% Magnitude part
subplot(2,1,1);
plot(f,abs(Y));
xlabel('Frequency (Hz)');
ylabel('Magnitude');
grid on;
% Phase part
subplot(2,1,2);
plot(f,unwrap(angle(Y))*180/pi); % in degrees (unwrapped)
xlabel('Frequency (Hz)');
ylabel('Phase (degrees)');
grid on;

If you zoom in on the peak in the magnitude spectrum, you will see that it is correctly at 377Hz.

I have put lots of comments in the code, but feel free to ask any questions.

 

[Saad Muftah]

I am really thankful to you, and thanks for the warm welcome, for the code it is look easy and that is because of your comments, i will check it and i ll ask you if there is any questions.

Many thanks,

 

[weetabixharry]

Sorry, I wrote this very quickly and made a minor error. Please note my correction (I have edited the code in my original post).

My error was that I defined t from zero to T_sine (including t = T_sine). My comments in the code only make sense now that I have defined:
t = (0:Nsamps-1)*T_sine/Nsamps;
T_samp = t(2)-t(1);

 

[Saad Muftah]

Thanks for the algorithm, but i have couple of questions:
1. why the magnitudes in the frequency domain are not the same as the signal it self, should it scaled.
2. and if you could help me in that, i want to build code extract magnitude of one frequency from signal that change with time, because the signal is transient and the harmonics change with time, so i want to see the behavior of one frequency with time.

 

[weetabixharry]

Thanks for the algorithm, but i have couple of questions:
1. why the magnitudes in the frequency domain are not the same as the signal it self, should it scaled.

Parseval's theorem states that the total energy in the time-domain signal, x(t), summed across all time t is equal to the total energy of its Fourier transform, Y(f), summed across all frequency components f. To check this in our discrete-time (Matlab) example, we can use the equation from this Wikipedia article:

https://en.wikipedia.org/wiki/Parseval's_theorem#Notation_used_in_engineering_and_physics

(the last equation before "See Also"). In Matlab, we write:

Energy_x = sum(abs(x).^2)
Energy_Y = (1/Nsamps)*sum(abs(Y).^2)

and find they are equal. Therefore, the scaling is correct.

2. and if you could help me in that, i want to build code extract magnitude of one frequency from signal that change with time, because the signal is transient and the harmonics change with time, so i want to see the behavior of one frequency with time.

Your question is not very clear. If you want to look at how frequency changes with time, a simple method would be to use a "sliding window" approach (e.g. if you have length(x) = 1000, you can compute Y1 = fft(x(1:16)); Y2 = fft(x(9:24)); Y3 = fft(x(17:32)); ... and so on). There are many better approaches, but this might get you started.

- - - Updated - - -

Another point I should probably add is that some people prefer to display zero Hertz (DC) at the middle of their spectrum. To do this, we can redefine Y and f as follows:

Y = fftshift(Y);
f = [-(ceil((Nsamps-1)/2):-1:1)*df 0 (1:floor((Nsamps-1)/2))*df];

For our simple sine example, this will reveal some nice symmetry (since purely real signals have conjugate symmetric spectra).

 

[Saad Muftah]

thanks for reply, for the second point in last reply; yes, it is something like that i want to absorb the magnitude of second harmonic of transient signal that decay with time, so my question how to make code for extracting the second harmonic in fourier spectrum. Another question if you don't mind how to get more knowledge in this field precisely, how to be more professional at this field, again many thanks,

 

[weetabixharry]

i want to absorb the magnitude of second harmonic of transient signal that decay with time, so my question how to make code for extracting the second harmonic in fourier spectrum.

I am still not sure exactly what you mean. When you say "absorb", do you mean remove? If you could provide more detail and an example (preferably with the data you want to process), I may be able to help.

Another question if you don't mind how to get more knowledge in this field precisely, how to be more professional at this field, again many thanks,

This is a large field, so I certainly wouldn't consider myself an expert. Personally, I just developed a bit of experience by working with examples in Matlab/C++ and searching the internet when I got stuck. However, I have a friend (who sat next to me during our PhDs) who is an expert and he highly recommends the following book for beginners:

"Who Is Fourier?: A Mathematical Adventure", Transnational College of LEX, 1995.

It is an unusual book, because it is written in a very easy-to-understand style, but actually covers the basics of Fourier quite thoroughly.

 

[Saad Muftah]

At the first thanks for the reference, and i want to explain my research task to you to full understand what iam looking for: it is a phenomena happened in transformer called inrush current. This current is transient current start very high and then decaying until reach steady state value, this phenomena is healthy phenomena because it is doesn't effect the transformer because it is temporary situation and then return normal. The digital protection relay could make false trip due to it, for that; i want to get algorithm that extract the second frequency" because the inrush current have a high value of second harmonic about 0.8 of the fundamental signal", so i want to extract that 2nd harmonic of this signal, and use it to notify the relay this is inrush current not fault or something harmful. One of those algorithm is by DFT or FFT or least square method.

Sorry i talk to long, but i want you to fully understand my task to enhance me.

Thanks;

 

[weetabixharry]

So, you want to compute the spectrum, then the frequency bin with the largest value tells you the principle frequency of the transient. Then, look at the frequency bin at double the frequency. If this is large enough (e.g., relative to the total energy in the spectrum), then everything is okay. If not, then trip.

 

[Saad Muftah]

Exactly, i don't know if you suggest particular algorithm or way to make such code. I will be apparatus if yo could you help me.
Many thanks,

 

[weetabixharry]

Well... do you have some real data (i.e. recordings of transients) you are working with? Or do you need to need to simulate some transients?

If you have some real data, then you need to load it first. Depending on the specific format, there are many ways of loading data in Matlab. For "raw" data, I use something like:

fid = fopen('filename.dat');
x = fread(fid,'float');

where "filename.dat" is the name of your data file, and the data is stored as single precision floating point numbers (i.e. 'float' ... but it can be doubles, or some kind of integers, etc.). If the data is formatted (not raw), you may need to use csvread() or dlmread().

If you need to simulate your own data, then I guess you could start with something like sine waves plus white noise (using the randn() function to produce additive white Gaussian noise). Or you can use a more detailed model if you want.

Then you can compute the spectrum, as I said before, with Y=fft(x). To find the maximum:

[max_val, max_loc] = max(abs(Y));
f_peak = f(max_loc);

So, f_peak is the principle frequency of the transient. Then, the next harmonic should be at frequency f(2*max_loc) with spectral magnitude abs(Y(2*max_loc)) (assuming your spectrum is big enough).

So, the basic information you have is: (1) Principle frequency (i.e. 1st harmonic), (2) Power of 1st harmonic (3) Frequency of 2nd harmonic (4) Power of 2nd harmonic. There are many ways you can use this information to do what you want. This is the fun part -- you can do some research and try different approaches. Maybe you don't even want to use the Fourier transform at all. There are many ways to approach this problem.

- - - Updated - - -

Also, as you move closer to designing a proper algorithm, and you have a larger amount of data to process, I would encourage you to consider replacing fft(x) with a better method for estimating Power Spectral Density. For example, if your transient exists during 33554432 data samples (i.e. 32*1024*1024 samples), but you only want to plot a spectrum with 8192 points, you could use Welch's sliding window method (e.g. in Matlab pwelch(x, gausswin(8192,5))).

If this is relevant to you, then I can try to provide more information. However, if your transients are very short, then this type of method may not be suitable.

 

[Saad Muftah]

Yes, the information you provided are very valuable for me, you put me on track. For the transient i am working on it's loaded from another program i don't need to generate any function, but the loaded signal is has a large number of data and as far as i know fft work with one period, and i need to do it for the whole signal, in other word i need to window it many times. And what you said about using other functions such as walsh and least square method is useful so much for me.
Can you help in this thing, and i would be appreciate.

Many Thanks my friend.

 

[weetabixharry]

First of all, I highly recommend these lecture notes on the theory of FFT-based Spectral Estimation: (click here). I think there is a lot of bad information on the internet about this, but those notes are very good. It may be a bit too theoretical for you, but it is still worth reading quickly.

In practice, we can start with a simple Matlab example. I have written some code quickly to hopefully show the basic idea of a sliding window method:

close all; clear all; clc;

% Total number of data samples in time domain
L = 1024*1024;

% Sliding window size in time domain (spectrum size in frequency domain)
Nfft = 512;

% Set the amount the windows overlap in the time domain
overlap = 0.5; % e.g. 50% overlap

% Generate some test data. Example: sine wave in Gaussian noise
x = sin(2*pi*(1/32)*(0:L-1)) + 10*randn(1,L);

% Plot a few sine periods (may look very noisy)
plot(x(1:128));
grid on;
axis tight;
title('Time-domain Input Data');


% Allocate storage for spectrum (i.e. Power Spectral Density estimate)
PSD = zeros(1,Nfft);
% Set window weights
w = gausswin(Nfft,5).'; % (e.g. a Gaussian window)
% Initialise sample index (i.e. where the window starts)
i = 1;
% Keep count of the number of FFTs that go into the PSD estimation
count_ffts = 0;
% Slide across all data
while i <= L - Nfft + 1

    % Select the next segment of data and multiply with window weights
    x_window = w.*x(i:i+Nfft-1);
    
    % Calculate power spectrum (unnormalised)
    P = abs(fft(x_window)).^2;
    
    % Add this power estimate to our total 
    PSD = PSD + P;
    count_ffts = count_ffts+1;

    % Increment sample index for next loop iteration
    i = i + round(Nfft*(1-overlap));
end

figure();
plot((0:Nfft-1),PSD);
title(['Unnormalised spectrum, averaged over ', num2str(count_ffts), 'spectra']);
grid on;
axis tight;

Note that I have not normalised the spectrum in any way, because this is often not necessary in practice. I don't think you will need to, unless your recording instruments were precisely calibrated and you need to measure specific quantities.

Of course, if you are looking for a transient, then one spectrum will not be enough; you will need the spectrum to update constantly. Therefore, you must decide how often you want the spectrum to be updated. This is an important trade-off. If the spectrum updates too slowly, then there may be a long delay between the transient occurring and being detected. However, if the spectrum updates too fast, it could have poor frequency resolution (Nfft too small) or be too noisy (not enough spectral averaging). You can experiment with different approaches.

 

[Saad Muftah]

I will start with that, and i will tell you if their is something stuck me.