1.760 views (last 30 days)

How can I find the frequency of this signal?

Thanks,

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

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.

Sign in to comment.

Image Analyst
on 25 Oct 2014

Image Analyst
on 26 Oct 2014

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

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

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:

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.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.