即以此功德,庄严佛净土。上报四重恩,下济三途苦。惟愿见闻者,悉发菩提心。在世富贵全,往生极乐国。
核心:czt可以在任何频段采样。假设同样的采样频率1kHz,同样的1000个点,如果FFT,你的分辨率是1Hz,但是如果你用czt,你可以在100和200Hz之间采1000个点,所以分辨率就变成了0.1Hz
clear; close all;
% Use czt to zoom in on a narrow-band section (100 to 150 Hz) of a filter's frequency response. First design the filter:
fs = 1e3;
Wc = 125; % cutoff frequency
h = fir1(30,125/(fs/2), rectwin(31)); % filter
% Establish frequency and CZT parameters.
% Restrict the frequency range of the CZT to the band between 100 and 150 Hz.
f1 = 100; f2 = 150; % in hertz
m = 1024; % FFT length, CZT length
% 2*pi*(f2-f1)/fs is the normalized angular frequency a.k.a. digital frequency, when f1 = 0, f2 = fs, this term has the maximum value 2*pi
w = exp(-j*2*pi*(f2-f1)/(m*fs)); % evenly divide 2*pi*(f2-f1)/fs into m parts
a = exp(j*2*pi*f1/fs); % starting angle 2*pi*f1/fs
z_spiral = a*(w.^-(0:m-1)), % the CZT spiral or "chirp" in the z-plane
figure; zplane([], z_spiral.');
z = czt(h,m,w,a);
% Compute the frequency response of the filter using fft and czt:
y = fft(h, m);
fn = (0:m-1)'/m;
fy = fs*fn;
fz = (f2-f1)*fn + f1;
subplot(211); plot(fy(1:m/2),abs(y(1:m/2))); axis([1 m/2 0 1.2]); xlabel('Hz'); ylabel('Magnitude'); title('Magnitude Response using FFT')
subplot(212); plot(fz,abs(z)); axis([f1 f2 0 1.2]); xlabel('Hz'); ylabel('Magnitude'); title('Magnitude Response using CZT ');
figure;
subplot(2,1,1);
plot(fy,abs(y));
axis([f1 f2 0 1.2]);
title('FFT');
subplot(2,1,2);
plot(fz,abs(z));
axis([f1 f2 0 1.2]);
title('CZT');
xlabel('Frequency (Hz)');
评论