How to get the exact amplitudes of an irregular wave?

I have time and displacement data of an irregular wave, and generated the amplitude spectrum by fft but the spectrum is very peaky it has lot of amplitudes (figure2). I need to extract the exact amplitude from that spectrum. For eg: given fig 1 has two peak amplitudes i need to get those amplitude at one time. If there was only one peak then it is easy to check by giving max command but for multiple peaks i dont know how to get that.

After fft i got the below figure, how to i get the exact amplitude values from the below figure?. Is there is any other way to extract the amplitudes?

The code i used for irregular sinwave to get the amplitude is,

T=time

X=Displacement

L = length(T);

nfft = 2^nextpow2(L);

f = Fs/2*linspace(0,1,nfft/2+1);

xfft = fft(X,nfft)/L;

xfftwant = 2*xfft(1:length(xfft)/2+1);

plot(f,abs(xfftwant ))

Answer 1

[Star Strider]

The second figure displays broadband noise. The easiest way to deal with that would be to denoise the original time-domain signal, either using the smoothdata function, the Wavelet Toolbox functions, or the Signal Processing Toolbox sgolayfilt function. The reason is that denoising the time-domain signal may be the easiest, since that allows interactive visual determination of the smoothed signal and the original signal. Then take the fft of the denoised signal to get the spectrum you want.

With broadband noise, frequency-selective filters will not work.

[Martin Thomas]

Hi,

I used the smoothdata , The graphs are smoothed actual wave form and amplitude spec. Is any chance to get the peakes of the amplitude spectrum. From the first question first graph there is only two peaks i need to get that peaks(both two top values) and the index at one time. Same can be applied for this too i think so. Please help me find a solution.

[Star Strider]

These seem to be the original data, not smoothed data. It is necessary to experiment with the smoothing to get the best results. Then do the analysis on the smoothed data.

[Star Strider]

I did some processing of your data, and I get significantly different results.

First, it is straightforward to read your data with textscan.

Second, the sampling intervals vary significantly, although they average out to produce a sampling frequency of about 128 Hz. After resampling them to a constant sampling interval (necessary for all signal processing including the fft), the data are much less noisy and relatively easy to interpret.

fidi1 = fopen('stillwaterT223H75.tom','rt');
D1C = textscan(fidi1, '%f%f%f%f', 'Delimiter','\t', 'CollectOutput',1);
D1 = cell2mat(D1C);
fclose(fidi1);

fidi2 = fopen('Te223Hs75.tom','rt');
D2C = textscan(fidi2, '%f%f%f%f', 'Delimiter','\t', 'CollectOutput',1);
D2 = cell2mat(D2C);
fclose(fidi2);

Fs=128;                 %sampling frequency
[D1r,t1r] = resample(D1(:,2:end), D1(:,1), Fs);                     % Resample To Constant Sampling Interval
[D2r,t2r] = resample(D2(:,2:end), D2(:,1), Fs);                     % Resample To Constant Sampling Interval

figure
plot(t1r, D1r)
grid
xlabel('Time')
ylabel('Amplitude')
title('stillwaterT223H75')

figure
plot(t2r, D2r)
grid
xlabel('Time')
ylabel('Amplitude')
title('Te223Hs75')

D2rs = D2r - mean(D1r);                                             % Subtract ‘mean(stillwaterT223H75)’ From ‘Te223Hs75’
D2rsm = mean(D2rs);

figure
plot(t2r, D2rs)
grid
title('Te223Hs75 Zeroed')

D2rss = smoothdata(D2rs, 'gaussian', 1E+3);                         % Choose Smoothing Method & Window

figure
plot(t2r, D2rss)
grid
title('Te223Hs75 Smoothed')

L = numel(t2r);                                                     % Data Length
Fn = Fs/2;                                                          % Nyquist Frequency
FTD2rs = fft(D2rss - mean(D2rss))/L;
Fv = linspace(0, 1, fix(L/2)+1)*Fn;                                 % Frequency Vector
Ix = 1:numel(Fv);                                                   % Index Vector

figure
nsp = size(FTD2rs,2);
for k = 1:nsp
    subplot(nsp,1,k)
    plot(Fv, abs(FTD2rs(Ix,k))*2)
    grid
    xlim([0 0.1])
    ylabel('Amplitude')
    title(sprintf('Column %d',k))
end
xlabel('Frequency (Hz)')

The Fourier transform results (after subtracting the means of each column to eliminiate the D-C offset and make the remaining peaks more apparent) are then:

I found your code somewhat difficult to follow, however I believe I captured the essence of it in my code.

Experiment with my code — especially with the smoothdata arguments — to get the result you want. I have no idea what I am looking at, or what the ‘correct’ result would be.

.

[Martin Thomas]

hey star, from the very first quest i have given the fft code i calculated the phase from xfft ie, phase=angle(xfft). But i am getting more value than the inuput value. why it happening so any solution to get exact count.

[Star Strider]

I have no idea what the problem is. The phase is bounded (unless unwrap is used with it) so that will be fixed. The amplitude of the Fourier transform in my code is slightly less that the amplitude of the time-domain ampliltude, simply beccause of the way the fft is calculated. My code is correct.

If the amplitude of the Fourier transform is greater than the amplitude of the time-domain signal, something is wrong in the code. I re-wrote your code because your code was so difficult to follow that I could barely make sense of it. It would be best to use my version.

[Martin Thomas]

that i understood, your code good but imtrying to follow the above one, beacuse that's my prof ask me to do. i have to sort it out my phase issue. no idea

[Star Strider]

The phase by definition is bounded by (or ).

I have no idea what the problem is, either, since I got good results.

Note that your data are not sampled uniformly, so that is the reason I used the resample function before I did my analysis. Non-uniformly sampled data will cause spurious Fourier transform results. Perhaps that is what you are seeing if you are not using my code.

I have no idea what you are doing, what your data are, or what the results should be. You need to ask your professor about these problems.

[Martin Thomas]

ya thats the reason, One thing i need to know is about my code with technicaly part

phase=angle(xfft(1:length(xfft)/2+1)); is this correct any technical error im getting a result with this but i need to know any technicaly worst idea or not.? because this will give the same lenth of xfftwant thats why. TO be frank is this a worst approch?

[Martin Thomas]

Yes that's the reason, One thing i need to know about my code from technical part.

I did the phase angle ike this

phase=angle(xfft(1:length(xfft)/2+1));

is this correct any technical error?

im getting a result but i need to know that is this a worst apporch or not . From this i will get the same lenth of amplitude and frequency.

[Martin Thomas]

hey star, the output was more than input length, now its correct, because i only need the + side.

[Star Strider]

Correct.

Only use the first half of the fft result (using the ‘Ix’ subscript range in my code) to plot the one-sided Fourier transform (actually Bode plot) displaying magnitude and phase.

[Martin Thomas]

from your code how can i seperate all the amplitudes in different columns. Your all in one abs(FTD2rs(Ix,k))*2 matrix. i need all the colums seperatly, i used (:,1) but its not working

[Star Strider]

See my previous Comment. I separated them in the subplot calls in the loop at the end.

This for loop separates them:

figure
nsp = size(FTD2rs,2);
for k = 1:nsp
    subplot(nsp,1,k)
    plot(Fv, abs(FTD2rs(Ix,k))*2)
    grid
    xlim([0 0.1])
    ylabel('Amplitude')
    title(sprintf('Column %d',k))
end

The ‘Ix’ vector selects the first half of the fft result, and ‘k’ selects the column.

Change the xlim arguments to see more of the plot, or different segments of it, or put a % in front of it temporarily to see the enitre spectrum.

[Martin Thomas]

hey star,

Time period =1/frequency.

Is that important to take the freq in terms of hz? I took in terms of radian/s.

[Star Strider]

It depends on what you want. The inverse of frequency is in uints of (time units)/(frequency units), so either seconds/cycle (original frequency Hz), or seconds/radian (original frequency radians/sec).

It should make no difference providing all the units are consistent throughout the code.

[Martin Thomas]

Hey star, one doubt before doing a numerical analysis in matlab with data, why the validation process is inmportant?

[Star Strider]

I don’t understand what ‘validation process’ refers to.

If my Answer helped you solve your problem, please Accept it!

.

[Martin Thomas]

hey one last doubt, fft result has more values than we need how to trim that values.? how to identify?

[Star Strider]

I am not certain what you want to do. I would not change the fft call to produce fewer points, although that is an option using the optional second argument to the fft function. Depending on what you want, you can just trim the end of the ‘FTD2rs’ matrix, or decimate it, for example as:

FTD2rsd = FTD2rs(1:10:end,:);

to keep only every tenth value. Change the ‘step’ value (10 here) to produce different results and different row-length ‘FTD2rsd’ variables. There are other possibilities as well, such as interpolating it to a new frequency vector, using linspace to define the interopolation vector and interp1 to do the actual interpolation. The choice is yours.