Finding the frequency value of a signal

How can I find the frequency of this signal?

Thanks,

[Yasemin G]

Hello Ramo,

Can you share your .m file with me please I am having the same problem.

[Image Analyst]

@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'

Answer 1

[Harry]

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']);

[sudi reddy]

I am new to matlab , when I use same code for spectrum , i didn't get any except that noisy signal.

Please let me know what are required to draw a spectrum

[Image Analyst]

pwelch() and plot() are required. We can't say anything more without your code and data.

[krn99]

What is f here, matlab throwing error saying 'undefined variable f'

[Image Analyst]

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]

From where NFFT = 1024 value came?

[Mahdi Saleh]

How could the center frequency be estimated in case you are dealing with modulated signal (QPSK, FSK, and others)? Is still this method valid?

Answer 2

[Image Analyst]

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.

[Ramo]

I smoothed the signal however couldn't find the peaks, i guess i got too much points, about 20000 points which form the signal. so the peaks are 3900 points.

[Image Analyst]

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]

>>

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]

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'

Answer 3

[Ramo]

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

[Harry]

Please upload your data file to this thread and I will modify my code to work with your data.

[Ramo]

attached

[Harry]

Ah okay, so the entire capture is just a few sine periods long. That's very short for an FFT-based method, since the resolution of the raw FFT is fsamp/Nfft (i.e. 10kHz resolution with your data).

I modified my code for your data and it gives an estimate of 40kHz. However, as I said, this is only accurate to the nearest 10kHz at best. Clearly, you could make an equally accurate estimate just by inspecting your graph (in your first post). Therefore, we should try a different approach.

One way is to launch Matlab's Curve Fitting Tool by typing:

sftool

Then, we type:

load uapp.txt
t = uapp(:,1);
x = uapp(:,2);

Then, inside the Curve Fitting Tool, we set "X data" and "Y data" as t and x. Finally, we select a "Sum of Sine" curve type, with "Number of terms" = 1. The result is:

The number we are interested in is the "b1 = 2.22e+05" on the left hand side. This means that the frequency estimate is:

2.22e+05/(2*pi) = 35332.3973664008 Hz

This should be significantly more accurate than the FFT-based method.

[Ramo]

I think I would need to use FFT as its fully automatic even though its not accurate, so if you could just modify your code to work with my signal. I need an automated method!

Thanks very much,

[Harry]

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
plot(F,t,x);
grid on;
xlabel('Time (secs)');
ylabel('Amplitude');
title(['Estimated frequency = ', num2str(F.b1/(2*pi)), ' Hz']);

This gives us a frequency estimate of about 35.33kHz:

[Ramo]

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]

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 ;