MATLAB Answers

Finding the frequency value of a signal

1.760 views (last 30 days)
Ramo on 25 Oct 2014
Commented: Image Analyst on 27 Mar 2020
How can I find the frequency of this signal?


Sign in to comment.

Accepted Answer

Harry on 26 Oct 2014
Edited: Harry on 26 Oct 2014
Whenever you're interested in frequency content of a signal, the Fast Fourier Transform is often an excellent tool to use (see help fft). More specifically, Matlab's PWELCH function will provide a Power Spectral Density estimate using Welch's method:
Here is an example of how to use it to estimate frequency:
close all; clear all; clc;
% Assume we capture 8192 samples at 1kHz sample rate
Nsamps = 8192;
fsamp = 1000;
Tsamp = 1/fsamp;
t = (0:Nsamps-1)*Tsamp;
% Assume the noisy signal is exactly 123Hz
fsig = 123;
signal = sin(2*pi*fsig*t);
noise = 1*randn(1,Nsamps);
x = signal + noise;
% Plot time-domain signal
plot(t, x);
ylabel('Amplitude'); xlabel('Time (secs)');
axis tight;
title('Noisy Input Signal');
% Choose FFT size and calculate spectrum
Nfft = 1024;
[Pxx,f] = pwelch(x,gausswin(Nfft),Nfft/2,Nfft,fsamp);
% Plot frequency spectrum
ylabel('PSD'); xlabel('Frequency (Hz)');
grid on;
% Get frequency estimate (spectral peak)
[~,loc] = max(Pxx);
title(['Frequency estimate = ',num2str(FREQ_ESTIMATE),' Hz']);


Show 2 older comments
krn99 on 6 Apr 2017
What is f here, matlab throwing error saying 'undefined variable f'
Image Analyst
Image Analyst on 6 Apr 2017
f is returned by pwelch() - see the documentation for that.
x is defined as "x = signal + noise;" so I don't know why it would say that, unless you tried to use it before you defined it.
Souarv De
Souarv De on 24 Sep 2019
From where NFFT = 1024 value came?

Sign in to comment.

More Answers (2)

Image Analyst
Image Analyst on 25 Oct 2014
I'd smooth it a bit with a 3rd order Savitzky-Golay filter, sgolayfilt() in the Signal Processing Toolbox, then I'd use findpeaks to get the period and 1/period is the frequency. Attached is a Savitzky-Golay filter demo.


Show 1 older comment
Image Analyst
Image Analyst on 26 Oct 2014
Please attach your signal data and code. You must not have smoothed it properly or used good parameters to findpeaks(). Try increasing the window length in sgolayfilt() until your signal doesn't have any noise.
Ramo on 26 Oct 2014
load uapp.txt
>> plot(uapp(:,1),uapp(:,2),'m-')
>> x = uapp(:,2);
>> smtlb = sgolayfilt(x,3,41);
plot(1:200000, x(1:200000))
axis([0 200000 -15000 15000])
title('normal graph')
% grid
axis([0 200000 -15000 15000])
title('smoothed graph')
% grid
Image Analyst
Image Analyst on 27 Mar 2020
Ramo, You chose a frame length, 41, that was far too small. It should be hundreds or thousands because you have 200 thousand data points (which is far too many to see on a single screen without zooming, by the way). Here is corrected code:
% Read in data.
data = dlmread('uapp.txt');
x = data(:, 1);
y = data(:, 2);
% Plot original noisy data.
subplot(1, 3, 2);
plot(x, y, '-');
title('Noisy Data', 'FontSize', 15);
grid on;
% Smooth the data
smoothedY = sgolayfilt(y, 4, 2001);
% Plot smoothed data.
subplot(1, 3, 3);
plot(x, smoothedY, 'r-', 'LineWidth', 2);
grid on;
title('Smoothed Data', 'FontSize', 15);
% Plot both original and smoothed data.
subplot(1, 3, 1);
plot(x, y, '-');
title('Noisy Data', 'FontSize', 15);
grid on;
hold on;
plot(x, smoothedY, 'r-', 'LineWidth', 2);
title('Both Noisy and Smoothed Data', 'FontSize', 15);
% Maximize the window.
g = gcf;
g.WindowState = 'maximized'

Sign in to comment.

Ramo on 26 Oct 2014
I got 20000 points (i.e. x and y) which i load from a txt file. then I plot them using plot(x,y).. what is the next step? Thanks


Show 4 older comments
Harry on 28 Oct 2014
I must reiterate that a basic FFT-based method is a very poor approach for such a short data capture (relative to the period of the sinewave), since it gives a very inaccurate result. You can even get a more accurate result just by looking at the graph and saying the period between the first peak and the second peak is about (40.2μs-12μs) = 28.2μs. Therefore, the frequency is about 1/28.2μs = 35.461kHz.
Instead, we can automate the curve-fitting method like this:
close all; clear all; clc;
% Load data
load uapp.txt
t = uapp(:,1);
x = uapp(:,2);
% Get sine fit
F = fit(t, x, 'sin1');
% Plot result
grid on;
xlabel('Time (secs)');
title(['Estimated frequency = ', num2str(F.b1/(2*pi)), ' Hz']);
This gives us a frequency estimate of about 35.33kHz:
Ramo on 30 Oct 2014
If you could explain the equation you have used please? and by the way I have tested your method with another signal it seems that it gave me the right frequency, however the sine wave didnt match at all.
Rodrigo Picos
Rodrigo Picos on 27 Mar 2020
Many years later.....
You should use 'sin3' or 'sin4', and check if you need to use the third or fourth component.
Hint: call F=fit( ) without the ending ;

Sign in to comment.

Sign in to answer this question.

Translated by