single sided amplitude fourier spectrum

Hello all, I am new user that trying to learn coding. I have an acceleration-time history with a time-step of 0.01 sec. The total duration is 14 sec. I would like to get a single-sided fourier spectrum of it. Is there any packed code for it? Any kind of help is appreciated. The attached file contains the data. Thank you.

Muhsin

Answer 1

[Star Strider]

See the documentation for fft (link).

The Code —

D = dlmread('Muhsin A00.txt', '\t', 1,0);
t = D(:,1);                                             % Time (s)
a = D(:,2);                                             % Acceleration (g)

L = length(t);
Ts = t(2)-t(1);                                         % Sampling Interval (sec)
Fs = 1/Ts;                                              % Sampling Frequency (Hz)
Fn = Fs/2;                                              % Nyquist Frequency (Hz)

FTa = fft(a)/L;                                         % Fourier Transform (Scaled)
Fv = linspace(0, 1, fix(L/2)+1)*Fn;                     % Frequency Vector
Iv = 1:length(Fv);                                      % Index Vector

figure(1)
plot(Fv, abs(FTa(Iv))*2)                                % One-Sided Amplitude Plot
xlabel('Frequency (Hz)')
ylabel('Amplitude (g)')
grid

EDIT — Changed ‘Ts’ calculation to accommodate ‘t(1)=0.01’ and ‘t(end)=0’.

[Muhsin]

Dear Star Strider, As far as I know single sided amplitude contains only the positive frequency. It discards the negative energy. The code you provide gives the negative frequencies with unacceptable magnitude of units(frequency is in the range of 0 ~ (-1.4) *10^-6 Hz). Thank you

[Star Strider]

My code displays only the positive frequencies. (See edited code.) My initial ‘Ts’ calculation assumes that your time data are monotonically increasing over the range, when in reality they are not. My edited code corrects for this problem.

I also provided you with one approach to correctly read your file.

The Plot —

EDIT — Added plot.

[Muhsin]

Dear Star;

The plot that your code gives is not compatible with the plot that I have. I do not know where the difference comes from.

[Muhsin]

Dear Star;

If you can look at the attached pic, you will see I got a plot that is close to plot I have for verification. I have used fft function that Matlab provides in the following link,

https://www.mathworks.com/help/matlab/ref/fft.html

Thank you

[Image Analyst]

Star's plot looks pretty much like your desired result. The differences might be only due to subsampling as it displayed it small on your different screens. So is it solved? If so, "Accept this answer".

[Muhsin]

 Star's plot has the same trend but different magnitudes. Since I do not know the difference I would not say it is like my desired result. Please take a look at my desired plot

[Star Strider]

My pleasure.

The difference is that your data are zero-padded for some reason. I used them as you posted them.

Try this (limiting the data to only the non-zero values of ‘t’):

D = dlmread('Muhsin A00.txt', '\t', 1,0);
t = D(:,1);                                             % Time (s)
a = D(:,2);                                             % Acceleration (g)

t_end = find(t == 0, 1, 'first');                       % Last Non-Zero ‘t’
idx = 1:t_end-1;                                        % Index For Non-Zero ‘t’

L = length(t(idx));
Ts = t(2)-t(1);                                         % Sampling Interval (sec)
Fs = 1/Ts;                                              % Sampling Frequency (Hz)
Fn = Fs/2;                                              % Nyquist Frequency (Hz)

FTa = fft(a(idx))/L;                                    % Fourier Transform (Scaled)
Fv = linspace(0, 1, fix(L/2)+1)*Fn;                     % Frequency Vector
Iv = 1:length(Fv);                                      % Index Vector

figure(1)
plot(Fv, abs(FTa(Iv))*2)                                % One-Sided Amplitude Plot
xlabel('Frequency (Hz)')
ylabel('Amplitude (g)')
grid

I could not determine if you wanted to retain the zeros or eliminate them. I eliminated them in this code. Note that it also affects the frequency vector. If you want to zero-pad your signal to increase the frequency resolution, use the second ‘NFFT’ argument to the fft function to do that in MATLAB:

D = dlmread('Muhsin A00.txt', '\t', 1,0);
t = D(:,1);                                             % Time (s)
a = D(:,2);                                             % Acceleration (g)

t_end = find(t == 0, 1, 'first');                       % Last Non-Zero ‘t’
idx = 1:t_end-1;                                        % Index For Non-Zero ‘t’

L = length(t(idx));
Ts = t(2)-t(1);                                         % Sampling Interval (sec)
Fs = 1/Ts;                                              % Sampling Frequency (Hz)
Fn = Fs/2;                                              % Nyquist Frequency (Hz)

NFFT = 2^nextpow2(t_end-1);
FTa = fft(a(idx), NFFT)/L;                              % Fourier Transform (Scaled)
Fv = linspace(0, 1, fix(L/2)+1)*Fn;                     % Frequency Vector
Iv = 1:length(Fv);                                      % Index Vector

figure(1)
plot(Fv, abs(FTa(Iv))*2)                                % One-Sided Amplitude Plot
xlabel('Frequency (Hz)')
ylabel('Amplitude (g)')
grid

This is the correct way to calculate and plot the Fourier transform of your signal. The Nyquist frequency is 50 Hz, the maximum resolvable frequency for your sampling frequency of 100 Hz.

[Muhsin]

Dear Star; Is it possible if want to calculate the single sided spectrum between specific lines. For Example,between 1504 - 2904 for the attached file. Thank you.

[Star Strider]

Simply read the file into your workspace (we’ve covered that already), then just choose the lines (rows) you want, and calculate the fft for those.

The problem is that when I read your file and then choose those lines:

D_main = dlmread('A00.txt', '\t', 1,0);
D = D_main(1504:2904, :);
t = D(:,1);                                             % Time (s)
a = D(:,2);                                             % Acceleration (g)

both ‘t’ and ‘a’ are completely zero. The Fourier transform of those data would tell you nothing.

[Muhsin]

Would you please check why the following script does not work? It was actually between the lines of 503-1903. I am doing something wrong. Thank you.

D_main = dlmread('A00.txt', '\t', 1,0);

D = D_main(503:1903, :);

t = D(:,1); % Time (s)

a = D(:,2); % Acceleration (g)

t_end = find(t == 0, 1, 'first'); % Last Non-Zero ‘t’

idx = 1:t_end-1; % Index For Non-Zero ‘t’

L = length(t(idx)); Ts = t(504)-t(503); % Sampling Interval (sec)

Fs = 1/Ts; % Sampling Frequency (Hz)

Fn = Fs/2; % Nyquist Frequency (Hz)

FTa = fft(a(idx))/L; % Fourier Transform (Scaled)

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

Iv = 1:length(Fv); % Index Vector figure(1)

plot(Fv, abs(FTa(Iv))*2) % One-Sided Amplitude Plot

xlabel('Frequency (Hz)')

ylabel('Amplitude (g)')

grid

[Star Strider]

Use this code instead:

D_main = dlmread('A00.txt', '\t', 1,0);
D = D_main(503:1903, :);
t = D(:,1);                                             % Time (s)
a = D(:,2);                                             % Acceleration (g)

L = length(t);
Ts = t(2)-t(1);                                         % Sampling Interval (sec)
Fs = 1/Ts;                                              % Sampling Frequency (Hz)
Fn = Fs/2;                                              % Nyquist Frequency (Hz)

NFFT = 2^nextpow2(L);
FTa = fft(a, NFFT)/L;                                   % Fourier Transform (Scaled)
Fv = linspace(0, 1, fix(L/2)+1)*Fn;                     % Frequency Vector
Iv = 1:length(Fv);                                      % Index Vector

figure(1)
plot(Fv, abs(FTa(Iv))*2)                                % One-Sided Amplitude Plot
xlabel('Frequency (Hz)')
ylabel('Amplitude (g)')
grid

[Muhsin]

I get the following errors with the latest code. Thank you

[Star Strider]

It worked for me. (I ran it and it worked in the complete code I posted.)

I just now ran it again and it again worked without problems. Since I can’t reproduce the error, I have no idea what the problem could be.

[Muhsin]

I found the problem. Just to let you know once I remove the headlines (first 3 lines), the code worked.

[Star Strider]

The file you originally posted only had one header line.

[Muhsin]

Dear Star; I would like to smooth the fourier spectrum. (not only single-sided spectrum but also the fourier spectrum). How do I do that? Can you help me please? Thank you.

[Star Strider]

You can window it using the spectrogram or fsst functions, or you can use Savitzky-Golay filtering (the sgolayfilt function), depending on the result you want.

Savitzky-Golay filtering might be easier, however getting ithe inverse Fourier transform from it would likely be impossible, even if you also filtered the phase information.

It’s probably best to just leave it as it is.

[Muhsin]

How about smoothdata function. Wouldn't that work well? Thank you

[Star Strider]

The smoothdata function would likely work. (It was introduced in R2017a, and since I do not know what version you have, I did not suggest it.)

You would use it on the absolute value of your single-sided Fourier transformed data.

Answer 2

[Image Analyst]

Isn't that what pwelch() does? (Signal Processing Toolbox required for pwelch).

[Image Analyst]

You might also like spectrogram() function.

Or the Signal Analyzer app on the Apps tab of the tool ribbon.