How is pspectrum() different from FFT() magnitude plotted in dB?

19 vues (au cours des 30 derniers jours)
Jay Vaidya
Jay Vaidya le 20 Déc 2020
Commenté : Mathieu NOE le 21 Déc 2020
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.

Réponse acceptée

Mathieu NOE
Mathieu NOE le 20 Déc 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
  2 commentaires
Jay Vaidya
Jay Vaidya le 20 Déc 2020
Modifié(e) : Jay Vaidya le 20 Déc 2020
I don't have more than 1 second of data but my sampling frequecny was quite high in exerpiements (~ 500 MSa/s). So I have ~ 50 ms of data but having ~1e7 data points. In any case, I don't understand if changing Fs is applicable to my case because I have limited data set, unlike the mathematical expression that you wrote above which can simulate the data for any given time.
Mathieu NOE
Mathieu NOE le 21 Déc 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;

Connectez-vous pour commenter.

Plus de réponses (0)

Tags

Produits


Version

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by