How to find the frequency and amplitude of an oscillating signal?
97 views (last 30 days)
Show older comments
Rohan Kokate
on 29 Nov 2023
Commented: Star Strider
on 1 Dec 2023
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!
0 Comments
Accepted Answer
Star Strider
on 29 Nov 2023
Edited: Star Strider
on 29 Nov 2023
One approach —
T1 = readtable('book1.xlsx', 'VariableNamingRule','preserve')
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
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_Magnitude = pks
Peak_Width = width*Fv(2)
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.
.
8 Comments
Star Strider
on 1 Dec 2023
As always, my pleasure!
Essentially, the more frequencies that exist in a signal, the more the total energy is distributed among them.
More Answers (1)
Paul
on 30 Nov 2023
Edited: Paul
on 30 Nov 2023
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.
See Also
Categories
Find more on Measurements and Feature Extraction in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!