How to do a fast Fourier transform (fft) in MATLAB to calculate the spectrum of data from a mat file 

You are welcome, I' glad that you found the video useful.

@@chihieule156 Thanks for your inquiry. What did you measure in time domain? Do you have two-dimensional data there are well, e.g. the (magnetic) force as a function of time and position?

And here is my answer: Thanks for asking. This is not a bug, it is a feature. "See also: en.wikipedia.org/wiki/Fourier_series#Common_forms The FFT calculates the double DC component a_0. The factor 2 is then considered as a divisor under a_0 in the formula. The DC component is then a_0/2. PS: The FFT becomes clearly better, if one omits the last point in the definition range and the sine function has then no jump when "putting the values one behind the other". See my version of the script:" x=linspace(0,2*pi,400); % remove last value x=x(1:end-1); y=10*sin(3*x)+2; figure plot(x,y) [f,X_f]=fourier(x,y,'sinus'); figure plot(f,abs(X_f))

Thanks for your question. By 'psd', you mean the 'power spectral density' (see: en.wikipedia.org/wiki/Spectral_density#Power_spectral_density), right? Maybe I could also record a video for this. To be honest, I have never really used the power spectral density before. Nevertheless, there are existing MATLAB functions readily available for this, see: de.mathworks.com/help/signal/ug/spectral-analysis.html

@@MathiasMagdowski Hopefully you had already worked with power spectral density (PSD). If you make a video on this it will certainly clarify. Thank you!

​@@giulianosalerno1203Maybe, but certainly not within the next weeks as I am busy with conferences, international teaching exchange programs and distinguished lecturer talks.

Thanks for your great question. The idea of specifying this additional number n as an optional argument to the fft function has something to do with how the fft works. In general, the Fourier transform is an integral transform. MATLAB is numerical software, so it deals with discrete sampled data. An integral over discrete sampled data becomes a sum, and an analytic Fourier transform becomes a discrete Fourier transform, or DFT. However, the sum within this DFT is quite expensive to compute, especially for large numbers of samples. This was especially a problem in earlier decades when computers were not that fast. Of course, mathematicians and engineers solved the problem and invented the Fast Fourier Transform, or FFT, which is essentially a DFT but uses a special and very fast algorithm to compute it. However, this fast algorithm only works when the number of samples is a power of 2, so you would normally choose n to be a power of 2. With computers 15 or 20 years ago, it would really make a difference if you did the DFT of 1000 samples or an FFT of 1024 samples. Today, you would not notice the difference in computation time, unless - you have really many samples (billions) - you have a slow computer (e.g. only a microcontroller) - you are trying to calculate many spectra in a short time. In this case, it might help to just cut some of the data or fill in additional zeros to get to the nearest power-of-2 number of samples. I hope this answers your question.

@@MathiasMagdowski Thank you for your quick response. As I understand from the Matlab documentation on fft(x,n): "If n is not specified, it defaults to the length of the input sequence". If this is true, then the length(x) is likely NOT a power of 2. It is not clear to me if the Matlab fft(x) function internally calculates the closest nextpow2 closest to the length(x) when n is not specified. Should this concern me or might it influence the results if n is not specified as a power of 2? Sometimes very different results are obtained based on the choice of n. For example, fft(x,512) may look different than fft(x,16384), or even if arbitrary values are used such as fft(x,12345). Is this because the dataset is truncated (or zeros added) thus results are expected to be different? Which is the best practice?: - not specifying n and let Matlab default to the length of the input sequence, or - specifying n as close to the nextpow2 as length(x) Thanks again!

@@louispasnik5976 Thank you for your further inquiry. If you leave n empty and the length of the input vector is not a power of 2, MATLAB will compute a general discrete Fourier transform that will give you the desired spectrum as a result, but the computation itself will not be "fast". In any case, MATLAB's fft function will NOT truncate or expand your input data to the closest length that is a power of 2. If you don't care about the speed or efficiency of your calculation, don't worry about n and just leave it as is. If you care about speed (because you have a huge vector, a very slow computer like a microcontroller, or if you are calculating many fft functions within a loop), try truncating or zero-padding your data to a length that is a power of 2. When you truncate data, you are of course losing data. If you zero-pad, you do not lose any information, but you still change your time function somehow, because it is no longer clear whether the function value is zero on purpose or just because of the padding. If you care about the speed or efficiency of your calculation, and if you have some control over the length of the input data (e.g., by setting a certain time step in a transient simulation, or by setting a certain sampling rate in a measurement), try to get something as close to a power of 2 as possible, so that you have to truncate or pad as few samples as possible. Of course, depending on x, you might get very different results for n set to 512, 16384, or something else, because many samples might be truncated and discarded, or added as zeros.

@@MathiasMagdowski Thank you, Professor!

@@louispasnik5976 You are welcome.

Thanks for your inquiry. The difference between the two modes 'pulse' and 'sinus' for the fourier.m function is also quite well documented within the source code of the function itself. The purpose of the 'pulse' is to calculate the Fourier transform of a single, non-repetitive pulse or an impulsive signal. In contrast, the idea of the 'sinus' mode is to calculate the Fourier transform of a continuous, repetitive signal like a pulse train. The 'pulse' mode corresponds to the common Fourier transform (en.wikipedia.org/wiki/Fourier_transform) whose result is a continous frequency spectrum or a density spectrum. If the time-domain function has for the example the unit V (volt), the density spectrum will have the unit V/Hz (volt per hertz). In contrast, the 'sinus' mode corresponds to the Fourier series (en.wikipedia.org/wiki/Fourier_series) whose result is a discrete frequency spectrum. Correspondingly, this spectrum will also have the unit V if the time function has the unit V. The main difference between the two modes inside the function 'fourier.m' is that a different scaling or renormalization of the amplitude is used.

If you obtain an error message stating that the argument data type is incorrect or the argument is missing, you should check if both if your argument variables t and s exist and what their data type is. You should be able to locate these variables in the MATLAB workspace and see their data type there. The data type should be a plain double-precision vector or array, not a cell array, table or something like this.

@@MathiasMagdowski thanks for your prompt reply. I figured it out. It had to do with way I loaded my data. I tweaked it a bit and everything was fine. Thank you very much 😊

Okay, very good. Here is a simple MATLAB code example that shows that for e.g. 2000 Hz sampling frequency, the highest frequency in the spectrum corresponds to 1000 Hz, just like one would expect due to the Nyquist-Shannon sampling theorem (en.wikipedia.org/wiki/Nyquist%E2%80%93Shannon_sampling_theorem): f_sample=2e3 f_signal=100 t_step=1/f_sample T_period=1/f_signal t=0:t_step:T_period-t_step x_t=5*cos(2*pi*f_signal*t) figure(1) plot(t,x_t) xlabel('time in s') ylabel('signal') [f,X_f]=fourier(t,x_t,'sinus') figure(2) plot(f,abs(X_f)) xlabel('frequency in Hz') ylabel('amplitude spectrum')

i can provide you more details with the generated data if needed.

@@sandeepvinnakota1552 Thanks for your inquiry and sorry for the delay in communication, but I was on vacation and on business travel. If you have a time-dependent function, e.g. the number of particles as a function of time, you can apply the Fourier transform and get the corresponding spectrum of this time function. This will allow you to analyse the frequency-dependent behavior of this process. For example, you may also observe dominant frequencies, if they exist. In contrast, I assume that it might be difficult to analyse, whether the particles are following a similar path, if this question is not already answered in your time-domain data. If you like, you may send me your well-described data set and your well-document MATLAB code to mathias.magdowski@ovgu.de and I might have a look.

Can I use linespace(start, end, no of points)?

If the timesteps are not equidistant, you need to interpolate them to equidistant time steps, e.g. using the interp1 function built into MATLAB, see: www.mathworks.com/help/matlab/ref/interp1.html You don't lose points or drastically lower the accuracy, if the equidistant (or uniform) time step of the interpolated time function is in the same range as the smallest time step of your non-uniform time steps. The linspace function is indeed useful to create equidistant time steps: t_uniform=linspace(0,t_end,N_steps); You can then interpolate your time function values: x_t_uniform=interp1(t_nonuniform,x_t_nonuniform,t_uniform); I hope this answers the question.

Thank you for your quick response. I tried to get the equidistant time step using the linspace() with start and end time but when I checked with plot(diff(time)) I didn’t get a line. I plotted the fft regardless but I am not sure the output (frequencies) is right. Would be great if you can have a look at my data.

@@ahmedbadawy9520 Sure, you may send me your data to my OVGU e-mail address.

Sorry for the delayed answer. I would suggest to run a for-loop over the channels and trials to calculate the corresponding spectrum separately. I'm not sure what you mean by the "spectrum of all the trials of each trials". Maybe you can explain this a little further. Is the average over all trials meant there or maybe the corresponding maximum? You may also send me your MATLAB code and some sample data via e-mail to: mathias.magdowski@ovgu.de

Thanks for the question. I don't have such video (yet), but I could record one for sure. Which part of the function is most unclear? There are also some comments in the fourier.m function. Which of them are difficult to understand?

hello,prof Thank you for your reply. In the function fourier.m, I have two questions: 1.the first one is " XX(1:N)=2/N*fft(x_t)", I can't understand it , for exemple, when I modify '2' to '2.5' in this sentence, the simulation of the image perfectly coincides with the correct image. 2.the second is " X_f(1:N/2+1)=XX(1:N/2+1)". I can't understand this sentence. Thank you very much. @@MathiasMagdowski

@@vivienne1592 Thanks for your reply. 1. If you change the factor 2 in the formula XX(1:N)=2/N*fft(x_t) to another value, e.g. 2.5, the shape of the spectrum will of course remain the same, because you change all the values of the array by a common factor. Still, the amplitudes given in the spectrum will not fit to the amplitude of the oscillations in the corresponding time function anymore. 2. The ordering of the frequencies in the output vector of the basic 'fft' function is kind of strange. The first entry in the spectrum obtained by the 'fft' function refers to the doubled constant part of the transient signal. It is followed by the positive frequencies in ascending order up the the highest frequency. Then, there are the negative frequencies in descending order. Negative frequencies only make sense for complex-valued time functions, which usually have no physical meaning in the engineering sciences. For real-valued time functions, the spectrum of the negative frequencies is the complex conjugate of the spectrum of the positive frequencies. Therefore, the negative frequencies can be omitted, as they are redundant and do not contain any relevant information in that case. Exactly this is done by the statement X_f(1:N/2+1)=XX(1:N/2+1). If I have time tomorrow, I might record a corresponding video with an explanation.

The requested video is now finished: "How to scale the result of the fft (fast Fourier transform) and how does fourier.m work in MATLAB?" ru-vid.com/video/%D0%B2%D0%B8%D0%B4%D0%B5%D0%BE-hAbsoXaZ144.html

Thank you very much! For 1, I can't understand, how to explain....my prof simulate the spectum with his script and I simulate the spectum with your script and my data.I known 2/N is a formula for normalisation... If I do not change the factor 2, the amplitude always different but I don't know if changing the factor from 2 to 2.5 is the correct or not, because in my simulation with the factor 2.5, my result plots and plots of my prof are completely consistent. @@MathiasMagdowski

Thanks for your question. Your task should not be too difficult to solve. If you want to get the amplitude and phase of the spectral component or fundamental frequency of your signal at a frequency of 60 Hz, your need to sample your time-domain signal with a sample rate of at least twice of this, so 120 Hz, but better with the 10-fold or even 100-fold of this frequency, so 600 Hz to 6000 kHz (corresponding to 6000 samples per second). You must also sample for at least one periodic time of this frequency of 60 Hz, so at least for 16.6 ms. Still, is is better to sample longer, up to 10 or even 100 periodic times, so 166.6 ms or even 1.666 s. Then you apply the FFT algorithm as explained in the video using the fourier.m function and just use the 60 Hz-component from the resulting frequency spectrum. I hope this helps. If you have further questions, please let me know. Best regards from Israel.

Thank you, another question , i want to calculate the current and voltage each sample to be able to detect when the disturbance occures for fault analysis, in another word i what to build for loop that calculate the impeadance from current and voltage each sample to monitor the change grahpically from the given rawdata

​@@sherifkhairy2879 Thanks for your further inquiry. Impedance is something that is only defined for (complex-valued) phasors of voltage and current, see en.wikipedia.org/wiki/Electrical_impedance Therefore, if you have a time-dependent voltage and time-dependent current function, you need to calculate their respective spectrum first. From the frequency-dependent voltage and current phasors, you can then calculate a frequency-dependent impedance. From this you might monitor some change in the behavior of your system for the detection of some occurence of disturbance or fault analysis.

Or this is realised using the introduced function automatically.

​@@hadinaderi7342 No, this is not done automatically by the provided fourier.m function. In general, for periodic functions (where the periodic time is well known and stable), it is good advice to always use a time range for the input of the FFT that is a multiple of the periodic time. You should keep in mind that the FFT assumes that the input signal always repeats or replicates itself consecutively. Therefore, the last time step of one time period should be omitted for the input so that the first time step of the next virtual replica can directly follow (without introducing a non-existing step, jump or sharp bend). In addition, one should consider that the time step of the time function corresponds to the highest frequency in the spectrum while the time range of the time function corresponds to the frequency step in the spectrum. I hope this helps.

Thanks for your question. Have you also used the fourier.m function from the ZIP file provided via the download link (cloud.ovgu.de/s/yrmZK2BCMW7t6KH) in the video description? This function takes two arguments, the time steps and the corresponding values of the time function. An optional third argument specifies the scaling of the spectrum. If you have another file or function present in your MATLAB workspace that is called fourier.m, this might lead to some confusion.

@@MathiasMagdowski yes, First I followed your tutorial with my own numbers and got the error. Than I downloaded your zip from under the video and got the same error. I can try the link you send me now later

@@fluminox13 Okay, feel free to come back to me at any time.

Thanks for your inquiry. After plotting, the y value or amplitude is still an absolute value in linear units, even if plotted on a logarithmic scale. If you want to convert this amplitude into decibels or dB, you need to do this conversion: en.wikipedia.org/wiki/Decibel#Root-power_(field)_quantities Does this answer your question?

@@MathiasMagdowski Thank you so much for your quick reply. Could you please check if my code is correct as follows? [f_CT,M_CT] = fourier(T,ID_CT_NG,'sinus'); plot(f_CT,20*log(abs(M_CT)),'blue'); Does this work? thank you so much!

@@xiadu3539 This works not really, as the decibel is defined using the decadic logarithm, i.e. the logarithm with the base of 10. The corresponding function is called log10 and not log in MATLAB.

Edit, Oh I am a beginner and didn't realise I had to put the source code files into the same folder as the written program. All sorted now, thank you very much! This has worked much better than the other fft tutorials, and also I am fft'ng a pulse signal so could use the pulse option, happy! :-)

@@Richard-ng2ch No problem, you are welcome.

Dear Tsehaye Negash, thanks for your further inquiry. There are some issues: 1. As I don't have your Microsoft Excel file "Addis_Ababa_wholeRS.xlsx", I cannot check your proposed MATLAB code. You may send it to me via mathias.magdowski@ovgu.de 2. Why do you not use the provided functions fourier.m and invfourier.m from cloud.ovgu.de/s/yrmZK2BCMW7t6KH, but try to use the more basic fft and ifft functions, where you have to take care about the frequency and time steps and the renormalization of the amplitudes yourself? Why make it difficult, if it can also be easy? 3. I don't really get your question? What shall the variable 'aAmplitude' mean or include? 4. Typically, one would like to plot amplitude vs. frequency, so the frequency on the x-axis or abscisssa and the amplitude on the y-axis or ordinate. Do you really want to plot the 'frequency vs amplitude'?

Each reference should have an author (Magdowski, Mathias), a title ("How to do a fast Fourier transform (fft) in MATLAB to calculate the spectrum of data from a mat file"), a year of publication (2021) and a way how to access it, e.g. a URL (ru-vid.com/video/%D0%B2%D0%B8%D0%B4%D0%B5%D0%BE-kUKqrd-lZBk.html or cloud.ovgu.de/s/yrmZK2BCMW7t6KH or ru-vid.com). I think this should be okay.

@@hariskhan485 Sure, you can do a Fourier transform for almost any time function, so also for market prices as function of time. The question is rather about the purpose of such analysis. As stock market prices behave rather random from my point of view, it will be difficult to perform an extrapolation or prediction for the future.

Dear Jason, thanks for your inquiry. You can find the corresponding fourier.m function in a ZIP file, which is linked in the video description: cloud.ovgu.de/s/yrmZK2BCMW7t6KH Does this help and answer your question?

Dear Mr. Dr. Magdowski, thank you very much for the fast respond it was helpful! With best regards Jason

@@jasonzimmermann9860 You are welcome.

Thanks for your question. As far as I understand your problem of particle image velocimetry (en.wikipedia.org/wiki/Particle_image_velocimetry), you have pictures with the positions of particles and want to calculate their velocity, which is the time derivative of the position. What do you need a Fourier transform or FFT for? Which frequency spectrum of which signal do you want to calculate?

@@MathiasMagdowski Yes, that is exactly what I meant by PIV. So, I have already processed the images and as a result I have a matrix of velocity values at different point in the picture. I want to use this matrix of velocity component to do a fft and plot it over time for any specific point in the image. I hope I was able to explain my problem properly.

@@mandarkadwekar6327 Thanks for the clarification. For me, it is still not clear what you would like to achieve. The result of the FFT will be a spectrum, i.e. the amplitudes and phases of sinusoidal signals with different frequencies. So the results of the FFT will be a function of frequency, not of time. If you want to do an FFT analysis and still plot something over time, you need to consider something like a short-time Fourier transform (en.wikipedia.org/wiki/Short-time_Fourier_transform). Nevertheless, there will always be a trade-off between frequency and time resolution.

@@MathiasMagdowski Hello thanks for your reply. Somehow I was confused while applying the code but now it is clear to me. I applied the source code now, but when I plot the function vs frequency graph, I am getting frequency for values that are higher than the maximum value in the original function matrix.

@@mandarkadwekar6327 Thanks for your further inquiry. I'm glad to hear that the MATLAB code is working for you. Nevertheless, I do not exactly understand your problem. What do you mean by "frequency for values that are higher than the maximum value in the original function matrix"? What is meant by "values", frequency values or amplitude/magnitude values? Again, as the results of the FFT or Fourier transform, you will get amplitudes and phases as a function of frequency, not the "frequency for values". What is further meant be the "original function matrix"? Why is it a matrix and not just a vector?

The 'fft' function (fast Fourier transform or discrete Fourier transform) is used within the 'fourier' function, which works like an API (application programming interface) or a wrapper to the underlying 'fft' function. The problems or issues with the basic 'fft' function are the following: 1. The 'fft' function only takes one argument as input, the values of the time function, and returns one output variable, the values of the spectrum. 2. Therefore, you have to take care about the corresponding time steps and frequency values yourself. The highest frequency corresponds to the shortest time, while the smallest frequency corresponds to the longest time. The 'fourier' function takes care of this. As a second input argument, you input corresponding time (equidistant) steps and get back the corresponding (equidistant) frequencies as a second output vector. 3. The ordering of the frequencies in the output vector of the basic 'fft' function is kind of strange. The first entry in the spectrum obtained by the 'fft' function refers to the doubled constant part of the transient signal. It is followed by the positive frequencies in ascending order up the the highest frequency. Then, there are the negative frequencies in descending order. Negative frequencies only make sense for complex-valued time functions, which usually have no physical meaning in the engineering sciences. For real-valued time functions, the spectrum of the negative frequencies is the complex conjugate of the spectrum of the positive frequencies. Therefore, the negative frequencies can be omitted, as they are redundant and do not contain any relevant information in that case. 4. The amplitude or magnitude of the spectrum calculated by the basic 'fft' function depends on the number of steps or vector length of the given time function. The wrapper provided by the 'fourier' function takes care of this and renormalizes the amplitude for two different cases, periodic functions (like sums of sine and cosine functions, rectangular pulse trains or clock signals) or aperiodic functions (like single pulses).

Thanks for the question. What do you exactly mean by a timestamp? Something like "2022-03-25 21:20:34.56"? If yes, you need to convert the timestamps into a vector of plain time steps with the unit second. Remember that the first time given in this time vector must be zero and that all other time steps must be equidistant. If they are not equidistant, you need to interpolate them to equidistant time steps, e.g. using the interp1 function.

I am not able to figure out how to convert it into a vector of seconds.. you think you could help with that?

@@kunalreddy434 I can try. Please send me your data file(s) and your current MATLAB code to my e-mail address mathias.magdowski@ovgu.de!

Sure! Thanks a lot

Any luck sir?