Finding the frequency value of a signal

Ramo (view profile)

on 25 Oct 2014
Latest activity Commented on by Souarv De

Souarv De (view profile)

on 24 Sep 2019 at 18:19
Accepted Answer by Harry

Harry (view profile)

How can I find the frequency of this signal?
Thanks, on 26 Oct 2014
Edited by Harry

Harry (view profile)

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

krn99 (view profile)

on 6 Apr 2017
What is f here, matlab throwing error saying 'undefined variable f'
Image Analyst

Image Analyst (view profile)

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 (view profile)

on 24 Sep 2019 at 18:19
From where NFFT = 1024 value came?

Answer by Image Analyst

Image Analyst (view profile)

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. Ramo

Ramo (view profile)

on 26 Oct 2014
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

Image Analyst (view profile)

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

Ramo (view profile)

on 26 Oct 2014
>>
>> 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

Ramo (view profile)

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

Ramo

Ramo (view profile)

on 28 Oct 2014
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

Harry (view profile)

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;
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

Ramo (view profile)

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. 