How to find the frequency and amplitude of an oscillating signal?

Hello,

I have a pressure signal which oscillates at a certain frequency and amplitude at steady-state. Is there a way to use the raw data (pressure and time) to find the frequency and amplitude of the oscillations. The raw data does have a little noise. I have attached an example image of the signal that I have along with an excel sheet which has the raw data. I believe FFT is the way to go but I have zero experience/background in using it.

Any help would be greatly appreciated.

Thank you in advance!

Answer 1

[Star Strider]

One approach —

T1 = readtable('book1.xlsx', 'VariableNamingRule','preserve')

T1 = 50×2 table
    Time [s]    P5 [psia]
    ________    _________

     3094.8      118.37  
     3097.1      116.52  
     3099.7      118.26  
     3102.2      120.05  
     3104.6      119.49  
     3107.1      116.54  
     3109.6      117.33  
     3112.1      119.37  
     3114.5      120.24  
       3117      117.57  
     3119.4      116.61  
     3121.9      118.29  
     3124.4      120.15  
     3126.8      119.49  
     3129.4      116.46  
     3131.8      117.45  

VN = T1.Properties.VariableNames;
t = T1{:,1};
s = T1{:,2};

figure
plot(t,s)
grid

checkTime = [mean(diff(t)) std(diff(t))]                                % Check Sampling Time Variation

checkTime = 1×2

2.4724 0.0631

Fs = 1/mean(diff(t));                                                   % Mean Sampling Frequency
[sr,tr] = resample(s,t,Fs);                                             % Resample To Constant Sampling Intervals (Required For Signal Ptocessing)
Fn = Fs/2;                                                              % Nyquist Frequency
L = size(T1,1);                                                         % Signal LEngth
NFFT = 2^nextpow2(L);                                                   % Use This Value To Make The 'fft' More Efficient
FTs = fft((s-mean(s)).*hann(L), NFFT)/L;                                % Discrete Windowed Fourier Transform
Fv = linspace(0, 1, NFFT/2+1)*Fn;                                       % Frequency Vector
Iv = 1:numel(Fv);                                                       % Index Vector

[pks,locs,width] = findpeaks(abs(FTs(Iv))*2, 'MinPEakProminence',0.25, 'WidthReference','halfheight');
Peak_Frequency = Fv(locs)

Peak_Frequency = 0.0885

Peak_Magnitude = pks

Peak_Magnitude = 0.8425

Peak_Width = width*Fv(2)

Peak_Width = 0.0168

idx = find(diff(sign(abs(FTs(Iv))*2-(pks/2))));
for k = 1:numel(idx)
    idxrng =  max(1,idx(k)-1) : min(idx(k)+1,Iv(end));
    hpf(k) = interp1(abs(FTs(idxrng))*2, Fv(idxrng), pks/2);
end

figure
plot(Fv, abs(FTs(Iv))*2)
hold on
plot(hpf, [1 1]*pks/2, '-k')
hold off
grid
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('Fourier Transform')
text(Peak_Frequency, Peak_Magnitude, sprintf('\\leftarrow Frequency = %.3f Hz\n    Magnitude = %.3f',Peak_Frequency, Peak_Magnitude), 'Horiz','left')
text(hpf(2), pks/2, sprintf('\\leftarrow FWHM = %.3f Hz', Peak_Width))

EDIT — (29 Nov 2023 at 18:28)

Added text calls to Forueir Transform plot.

.

[Rohan Kokate]

That makes sense! Thank you.

[Star Strider]

As always, my pleasure!

[Rohan Kokate]

I still am a little confused about the interpretation of the results from FFT. What is the difference between the magnitude of the FFT (= 0.842 units) and the amplitude of the oscillations (which I calculated manually from the data = 3.8 psia). How are these two parameters related? Essentially, my question is shouldn't these two values be the same?

[Star Strider]

The fft amplitude roughly corresponds to the peak-to-peak amplitude of the oscillations, with reservations for the fact that the peaks do not all have the same amplitudes, and the total energy of the waveform is distributed over about 0 to 0.2 Hz, and the peak itself has a full-width-half-maximum value of 0.017 Hz, meaning that the energy of the oscillation is distributed over that range. (That implies that even with the resampled signal, the oscillation frequency varies over a reasonable range.)

For a pure sine or cosine curve, the fft amplitude is almost exactly the peak-to-peak amplitude, and the width of the peak approaches zero. For real-world signals, it is usually different, and sometimes significantly so.

T1 = readtable('book1.xlsx', 'VariableNamingRule','preserve')

T1 = 50×2 table
    Time [s]    P5 [psia]
    ________    _________

     3094.8      118.37  
     3097.1      116.52  
     3099.7      118.26  
     3102.2      120.05  
     3104.6      119.49  
     3107.1      116.54  
     3109.6      117.33  
     3112.1      119.37  
     3114.5      120.24  
       3117      117.57  
     3119.4      116.61  
     3121.9      118.29  
     3124.4      120.15  
     3126.8      119.49  
     3129.4      116.46  
     3131.8      117.45  

VN = T1.Properties.VariableNames;
t = T1{:,1};
s = T1{:,2};

Fs = 1/mean(diff(t));                                                   % Mean Sampling Frequency
[sr,tr] = resample(s,t,Fs);                                             % Resample To Constant Sampling Intervals (Required For Signal Ptocessing)

[vy,pk] = bounds(sr);

figure
plot(tr,sr, 'DisplayName','Signal')
yline(pk, '--g', 'DisplayName','Maximum Peak Amplitude')
yline(vy, '--r', 'DisplayName','Minimum Valley Amplitude')
grid
legend('Location','best')
xlabel(VN{1})
ylabel(VN{2})
pos = get(gca,'Position');
xapf = @(x,pos,xl) pos(3)*(x-min(xl))/diff(xl)+pos(1);                % 'x' Annotation Position Function
yapf = @(y,pos,yl) pos(4)*(y-min(yl))/diff(yl)+pos(2);                % 'y' Annotation Position Function
annotation('doublearrow', xapf([1 1]*3090, pos, xlim), yapf([pk vy], pos,ylim))
text(3091, mean([pk vy]), sprintf('Peak-to-Peak Maximum = %.3f',pk-vy), 'Horiz','left', 'Vert','middle')

.

[Rohan Kokate]

I see. So essentially, the FFT magnitude does not really have any physical signficance for a real-world signal. The FFT magnitude only describes the output qualitatively. Is that a correct assesment?

[Star Strider]

Not really. The Fourier transform calculates the energy in a signal by frequency, and if the frequency varies, that will produce a wide peak with the energy of the signal dispersed over those frequencies. Add to that varying amplitudes, and the Fourier transform will do its best to describe the signal taking all that into account.

I devised a signal to demonstrate that and did a similar analysis on it here. Note that ‘s’ has an amplitude of approximately ±1 throughout its range, however the frequency varies significantly (50 Hz, with a centre frequency of 150 Hz), and the Fourier transform has an amplitude that is less tha 1% of the original signal amplitude. Some of this is ‘spectral leakage’, however most of it is due to the fact that the signal energy is spread over ±25 Hz from the centre frequency. That dissipates the energy significantly.

Example —

Fs = 1000;
Ls = 120;
t = linspace(0, Fs*Ls, Fs*Ls+1).'/Fs;
s = sum(sin(2*pi*t.*(100+50*t/max(t))),2);

figure
plot(t, s)
grid
xlabel('t')
ylabel('Amplitude')
title('Whole Signal')

figure
plot(t, s)
grid
xlabel('t')
ylabel('Amplitude')
title('Zoomed Segment')
xlim([-0.025 0.025]+60)

Fn = Fs/2;                                                              % Nyquist Frequency
L = numel(t);                                                           % Signal LEngth
NFFT = 2^nextpow2(L);                                                   % Use This Value To Make The 'fft' More Efficient
FTs = fft((s-mean(s)).*hann(L), NFFT)/L;                                % Discrete Windowed Fourier Transform
Fv = linspace(0, 1, NFFT/2+1)*Fn;                                       % Frequency Vector
Iv = 1:numel(Fv);                                                       % Index Vector

[pks,locs,width] = findpeaks(abs(FTs(Iv))*2, 'MinPeakProminence',0.005, 'WidthReference','halfheight');
Peak_Frequency = Fv(locs);
Peak_Magnitude = pks;
Peak_Width = width*Fv(2);

idx = find(diff(sign(abs(FTs(Iv))*2-(max(pks)/2))));
for k = 1:numel(idx)
    idxrng =  max(1,idx(k)-1) : min(idx(k)+1,Iv(end));
    hpf(k) = interp1(abs(FTs(idxrng))*2, Fv(idxrng), max(pks)/2);
end

figure
plot(Fv, abs(FTs(Iv))*2)
hold on
plot([min(hpf) max(hpf)], [1 1]*max(pks)/2, '-k')
hold off
grid
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('Fourier Transform')
text(Peak_Frequency, Peak_Magnitude, sprintf('\\leftarrow Frequency = %.3f Hz\n    Magnitude = %.3f',Peak_Frequency, Peak_Magnitude), 'Horiz','left')
text(hpf(2), pks/2, sprintf('\\leftarrow FWHM = %.3f Hz', Peak_Width))

.

[Rohan Kokate]

Oh, I get it now. I will try to read a bit more to get more clarity. Thanks for all the help!

[Star Strider]

As always, my pleasure!

Essentially, the more frequencies that exist in a signal, the more the total energy is distributed among them.

Answer 2

[Paul]

Hi Rohan,

For this signal, the output of fft without windowing gives a very good estimate of the amplitude and frequency of the signal.

Read and plot the data.

T = readtable('book1.xlsx', 'VariableNamingRule','preserve');
t = T{:,1};
p5 = T{:,2};
figure
plot(t,p5),grid

Get a rough estimate of its amplitude and frequency

[highpeaks,locs] = findpeaks(p5);
lowpeaks = -findpeaks(-p5);
amplitudeEstimate = (mean(highpeaks)-mean(lowpeaks))/2;
frequencyEstimate = 1./mean(diff(t(locs)));

The data is not equally spaced, but we'll assume it is for a first cut and use the mean dt to define a sampling frequency. Not necessarily a recommended approach, but may be good enough for a first cut because the deviation is not too large.

plot(1./diff(t))

Fs = 1/mean(diff(t));

Compute the FFT

P5 = fft(p5-mean(p5))/numel(p5);
f = (0:numel(P5)-1)/numel(P5)*Fs;

If you want to use a window, the scaling has to be modified

whann = hann(numel(p5));
P5hann = fft((p5-mean(p5)).*whann)/sum(whann);

Plot the results

plot(f,2*abs(P5),'-o',f,2*abs(P5hann),'-o'),grid
yline(amplitudeEstimate,'m'),xline(frequencyEstimate,'m')
xlim([0 Fs/2]),xlabel('Frequency (Hz)'),ylabel('Amplitude')
legend('Raw','Windowed')

You can probably sharpen this analysis by being more careful in dealing with the non-equally spaced samples, and maybe zero padding the fft's as well.

[Rohan Kokate]

Thanks Paul!