21 views (last 30 days)

Show older comments

I am having difficulty in finding a high-resolution frequency spectrum even though I have a pretty large number of time domain data samples (1999999 data points). The fundamental frequency peak in FFT is with a resolution of 50 Hz when I use fft(). I had also put this post here with the screenshots, code, and data.

I am thinking of using pspectrum() or zoomFFT(). The latter seems to be more relevant to me. But after creating the object of dsp.zoomFFT(), I am not sure how to proceed because z = zfft(data(:,2)) gives me error of

Fs = 40e3;

CF = 4e3;

BW = 1e3;

D = Fs/BW;

fftlen = 64;

L = D * fftlen;

zfft = dsp.ZoomFFT(D,CF,Fs,'FFTLength',fftlen);

z = zfft(data(:,2));

The number of rows in the input signal must be a multiple of 48 (the decimation factor).

Error in plot_fft_v8 (line 124)

Also, I am thnking that probably pspectrum() can help because

pspectrum(xTable,'FrequencyResolution',10)

Error using pspectrum>validateFrequencyResolution (line 1005)

'FrequencyResolution' value is not achievable for the specified parameters. 'FrequencyResolution' value must

be within [8.1702e-07, 1.634] (*pi radians/sample).

Error in pspectrum>parseAndValidateInputs (line 811)

validateFrequencyResolution(opts.FrequencyResolution, opts.Leakage,...

Error in pspectrum (line 246)

opts = parseAndValidateInputs(x,varargin);

can set my frequency resolution.

Mathieu NOE
on 20 Dec 2020

hello

if your data lenght is very long, there is no reason why you could not getter better resolution with a standard fft

remember that freq resolution df = Fs/NFFT

so if like in this demo I choose NFFT = Fs = 40e3, my resolution is 1 Hz, still you could even go beyong as you have more than 1 second of data (nb of samples is > Fs)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% load signal

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% dummy data

Fs = 40e3;

samples = 1999999;

t = (0:samples-1)'*1/Fs;

signal = cos(2*pi*50*t)+cos(2*pi*100*t)+1e-3*rand(samples,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% FFT parameters

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

NFFT = 40e3; %

Overlap = 0.75;

w = hanning(NFFT); % Hanning window / Use the HANN function to get a Hanning window which has the first and last zero-weighted samples.

%% notch filter section %%%%%%

% H(s) = (s^2 + 1) / (s^2 + s/Q + 1)

fc = 50; % notch freq

wc = 2*pi*fc;

Q = 5; % adjust Q factor for wider (low Q) / narrower (high Q) notch

% at f = fc the filter has gain = 0

w0 = 2*pi*fc/Fs;

alpha = sin(w0)/(2*Q);

b0 = 1;

b1 = -2*cos(w0);

b2 = 1;

a0 = 1 + alpha;

a1 = -2*cos(w0);

a2 = 1 - alpha;

% analog notch (for info)

num1=[1/wc^2 0 1];

den1=[1/wc^2 1/(wc*Q) 1];

% digital notch (for info)

num1z=[b0 b1 b2];

den1z=[a0 a1 a2];

freq = linspace(fc-1,fc+1,200);

[g1,p1] = bode(num1,den1,2*pi*freq);

g1db = 20*log10(g1);

[gd1,pd1] = dbode(num1z,den1z,1/Fs,2*pi*freq);

gd1db = 20*log10(gd1);

figure(1);

plot(freq,g1db,'b',freq,gd1db,'+r');

title(' Notch: H(s) = (s^2 + 1) / (s^2 + s/Q + 1)');

legend('analog','digital');

xlabel('Fréquence (Hz)');

ylabel(' dB')

% now let's filter the signal

signal_filtered = filtfilt(num1z,den1z,signal);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% display : averaged FFT spectrum before / after notch filter

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[freq,fft_spectrum] = myfft_peak(signal, Fs, NFFT, Overlap);

sensor_spectrum_dB = 20*log10(fft_spectrum);% convert to dB scale (ref = 1)

[freq,fft_spectrum_filtered] = myfft_peak(signal_filtered, Fs, NFFT, Overlap);

sensor_spectrum_filtered_dB = 20*log10(fft_spectrum_filtered);% convert to dB scale (ref = 1)

figure(2),semilogx(freq,sensor_spectrum_dB,'b',freq,sensor_spectrum_filtered_dB,'r');grid

title(['Averaged FFT Spectrum / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);

legend('before notch filter','after notch filter');

xlabel('Frequency (Hz)');ylabel(' dB')

function [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)

%FFT peak spectrum of signal (example sinus amplitude 1 = 0 dB after fft).

% signal - input signal,

% Fs - Sampling frequency (Hz).

% nfft - FFT window size

% Overlap - buffer overlap % (between 0 and 0.95)

samples = length(signal);

% fill signal with zeros if its length is lower than nfft

if samples<nfft

s_tmp = zeros(nfft,1);

s_tmp((1:samples)) = signal;

signal = s_tmp;

end

% window : hanning

window = hanning(nfft);

window = window(:);

% compute fft with overlap

offset = fix((1-Overlap)*nfft);

spectnum = 1+ fix((samples-nfft)/offset); % Number of windows

% % for info is equivalent to :

% noverlap = Overlap*nfft;

% spectnum = fix((samples-noverlap)/(nfft-noverlap)); % Number of windows

% main loop

fft_spectrum = 0;

for i=1:spectnum

start = (i-1)*offset;

sw = signal((1+start):(start+nfft)).*window;

fft_spectrum = fft_spectrum + (abs(fft(sw))*4/nfft); % X=fft(x.*hanning(N))*4/N; % hanning only

end

fft_spectrum = fft_spectrum/spectnum; % to do linear averaging scaling

% one sidded fft spectrum % Select first half

if rem(nfft,2) % nfft odd

select = (1:(nfft+1)/2)';

else

select = (1:nfft/2+1)';

end

fft_spectrum = fft_spectrum(select);

freq_vector = (select - 1)*Fs/nfft;

freq_vector = freq_vector(:);

fft_spectrum = fft_spectrum(:);

end

Mathieu NOE
on 21 Dec 2020

hello

from the attached picture, it seems that you are looking for a peak in the 4 kHz

your sampling at 500MS/s is overkill , regarding sampling (Shannon), 20 kS/s would suffice

that should allow you to increase the acquisition time by the same amount , so you could quite significantly improve your frequency resolution without having a huge file size;

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

Start Hunting!