Finding the frequency value of a signal

781 views (last 30 days)
How can I find the frequency of this signal?
Thanks,
  2 Comments
Yasemin G
Yasemin G on 26 Oct 2021
Hello Ramo,
Can you share your .m file with me please I am having the same problem.
Image Analyst
Image Analyst on 26 Oct 2021
@Yasemin G try this:
% Initialization steps.
clc; % Clear the command window.
close all; % Close all figures (except those of imtool.)
imtool close all; % Close all imtool figures if you have the Image Processing Toolbox.
clear; % Erase all existing variables. Or clearvars if you want.
workspace; % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 18;
%==================================================================================
% Image Analyst's code below:
period = 0.57e5
originalFrequency = 1/period
x = linspace(0, 2e5, 2000);
signalAmplitude = 0.85;
perfectSineWave = signalAmplitude * sin(2 * pi * (x - 0.08e5) / period);
subplot(2,1,1);
plot(x, perfectSineWave, 'b-');
grid on;
noiseAmplitude = 0.05;
yNoisy = perfectSineWave + noiseAmplitude * (2 * rand(1, length(x)) - 1);
hold on;
darkGreen = [0, 0.5, 0];
plot(x, yNoisy, '-', 'Color', darkGreen)
yline(0, 'Color', 'b', 'LineWidth', 2)
xlabel('Time or x', 'FontSize',fontSize);
ylabel('Signal', 'FontSize',fontSize);
title('Original Signal', 'FontSize',fontSize);
%==================================================================================
% Harry's code below:
% Assume we capture 8192 samples at 1kHz sample rate
Nsamps = 8192;
fsamp = 1000;
Tsamp = 1/fsamp;
t = (0:Nsamps-1)*Tsamp;
% Choose FFT size and calculate spectrum
Nfft = 1024;
[Pxx,f] = pwelch(yNoisy, gausswin(Nfft), Nfft/2, Nfft,fsamp);
% Plot frequency spectrum
subplot(2,1,2);
plot(f, Pxx, 'b-', 'LineWidth', 2);
ylabel('PSD', 'FontSize',fontSize);
xlabel('Frequency (Hz)', 'FontSize',fontSize);
grid on;
% Get frequency estimate (spectral peak)
[~,loc] = max(Pxx);
FREQ_ESTIMATE = f(loc)
caption = sprintf('Frequency estimate = %f Hz', FREQ_ESTIMATE);
title(caption, 'FontSize',fontSize);
g = gcf;
g.WindowState = 'maximized'

Sign in to comment.

Accepted Answer

Harry
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:
[Pxx,F] = pwelch(X,WINDOW,NOVERLAP,NFFT,Fs)
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
subplot(2,1,1);
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
subplot(2,1,2);
plot(f,Pxx);
ylabel('PSD'); xlabel('Frequency (Hz)');
grid on;
% Get frequency estimate (spectral peak)
[~,loc] = max(Pxx);
FREQ_ESTIMATE = f(loc)
title(['Frequency estimate = ',num2str(FREQ_ESTIMATE),' Hz']);
  5 Comments
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.
  4 Comments
Ramo
Ramo on 26 Oct 2014
>>
load uapp.txt
>> plot(uapp(:,1),uapp(:,2),'m-')
>> x = uapp(:,2);
>> smtlb = sgolayfilt(x,3,41);
subplot(2,1,1)
plot(1:200000, x(1:200000))
axis([0 200000 -15000 15000])
title('normal graph')
% grid
subplot(2,1,2)
plot(1:200000,smtlb(1:200000))
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
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
  7 Comments
Ramo
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.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!