How do I do the fourier transform of a data array?

I have been trying to do the fourier transform of a reference signal and of a row of signals, but the results I am getting are not the way they should be. What am I doing wrong?

I have got a two data arrays, the first is a reference signal and the second are 240 individual main signals (obtained over a period of 60 seconds). These signal are THz signals and their Fourier Transform is expected to take a certain shape, as shown in the the figure below (left shows the signal and right shows the fourier of the signal):

Fig.1 - Stereotypical FFT of a THz signal.

Each signal (the reference and the 240 main data signals) is composed of 1496 points.

The reference is : 1 x 1496;

The main signal is : 1496 x 240 (I was too lazy to make a transpose of this, but it can be easily done).

To get the fourier transform of said signal I have written the following:

ref_ft = fft(ref_clean);        %ref_clean is the reference signal 1 x 1496;
for p=1:240
    data_ft(:,p)=fft(data_clean(:,p)); %data clean are the main 240 signals, 1496 x 240 ;
end 
plot(freq, ref_clean);   %plotting them against frequency
plot(freq, data_clean)

Fig. 2 - My FFT of the reference signal

This, in no shape of form looks like examples shown in Fig. 1. What could I be possibly doing wrong? Is it potentially the shape of my data? I don't think it is because my reference is a simple shape and I haven't had to put it through a for loop.

Any help would be deeply appreciated.

I have attached the data files in case anyone finds something weird about them.

Answer 1

[Star Strider]

Try this —

% LD1 = load('ref_clean.mat');

% ref_clean = LD1.ref_clean;

LD2 = load('data_clean.mat');

data_clean = LD2.data_clean;

figure

waterfall(data_clean)

grid on

xlabel('rows')

ylabel('cols')

title('Signal Matrix')

Fs = 1; % Default Sampling Frequency (Actual Value Unstated)

Fn = Fs/2; % Nyquist Frequency

L = size(data_clean,1); % Signal Length

NFFT = 2^nextpow2(L); % For Efficiency

FTd_c = fft(data_clean,NFFT)/L;

Fv = linspace(0, 1, NFFT/2+1)*Fn; % Frequency Vector

Iv = 1:numel(Fv); % Index Vector

figure

for k = 1:fix(size(FTd_c,2)/20)

idx = 12*(k-1)+1;

subplot(6,2,k)

plot(Fv, abs(FTd_c(Iv,idx))*2)

grid

xlabel('Frequency (Units)')

ylabel('Amplitude')

title(sprintf('Row %3d',idx))

end

Fh = gcf;

pos = Fh.Position;

Fh.Position = pos + [-100 -500 200 500]; % Expand To Show Detail

I have no idea what ‘ref_clean’ is, so I didn’t use it.

All the signals appear to be essentially the same, so their Fourier transforms are also essentially the same.

Supply the appropriate value of ‘Fs’ to get the correct frequency vector.

.

[Goncalo Costa]

Thank you so much for your help. I have been trying to understand what you have written about, and have ran the code, and it works quite well, but I am still not certain as to how, and the reason for it to be written the way it is. It would be of great help if you could answer some of my questions. Thanks in advance.

• Why is the FFT divided by the length of the signal? I see that this has an impact on the amplitude of the signal, but I don't understand the reason.

FTd_c = fft(data_clean,NFFT)/L;

• Where does this formula for the frequency vector come from? I don't understand why it is multiplied by the Nyquist frequency, why the original limits are 0 to 1, or why this range is divided by NFFT/2+1 .... I am sorry, I know this is quite a lot.

Fv = linspace(0, 1, NFFT/2+1)*Fn;

Why are these the parameters chosen for the for loop? Is this simply because you were just wanting to show that it worked for the different rows of data? Because, otherwise, why not just put

t = 0:60/1496:60; % time-array
deltaT= t(2)- t(1);
MaxFreq=1/deltaT;
DeltaFreq = MaxFreq/(size(data_clean,1)-1);

freq=0:DeltaFreq:MaxFreq;
for k =1:width(data_clean)
    plot(freq, fft(data_clean(:,k), NFFT)/L)
end

I am sorry byt I don't understand the plot(Fv, abs(FTd_c(Iv,idx))*2) line very well.

Thank you very much for your help in advance.

[Star Strider]

My pleasure!

The fft result is divided by the length of the signal in order to normalise its amplitude. This is discussed more fully in the fft documentation.

That version of the frequency vector first appeared in the R2015b fft documentation, and I have used it ever since. (It is no longer available online. There are different expressions for it in different releases, for example beginning at least in R2017a it is ‘Fs*(0:(L/2))/L’ where ‘L’ is the length of the fft result [not necessarily the length of the signal since it has to accommodate the use of ‘NFFT’ if that is chosen], and ‘Fs’ is the sampling frequency.)

The loop to display some of the fft results is simply arbitrary. I chose a representative sample of 12 of them (there are 240 total) to display, and then stretched the resulting figure to make them easier to see.

This line:

plot(Fv, abs(FTd_c(Iv,idx))*2)

plots the frequency vector ‘Fv’ and the column of ‘FTd_c’ given by ‘idx’, and the index vector ‘Iv’ (to plot a one-sided fft result). The absolute value calculates the amplitude of the complex fft, and since the fft result is ‘symmetrical’, dividing the energy of the signal between the ‘positive’ and ‘negative’ frequencies (more easily appreciated when using fftshift), multiplied by 2 to compensate for displaying only the ‘positive’.frequencies in order to approximate tha amplitude of the time domain signal.

.

[Goncalo Costa]

Thank you very much for your help. I think I understand it now (at least most of it).

Each 240 measurements were taken over a period of 60 secs, so I changed Fs to be 4. I think this was the correct thing to do.....

I have also made alterations in the plot to plot every single Fourier transform.

For anyone else having doubts regarding the Frequency vector and the Index Vector, as well as the NFFT there are some links below that might help:

https://uk.mathworks.com/matlabcentral/answers/532448-help-me-undersandand-frequency-vector-and-index-vector-in-fourier-transform

https://uk.mathworks.com/matlabcentral/answers/324316-what-is-nfft-and-where-it-is-used-in-matlab-code

[Star Strider]

As always, my pleasure!