标准DFT/FFT的局限性

数字信号的光谱分辨率由它的长度决定(因此是有限的)。我们将用一个简单的例子来说明这个事实。考虑一个由两个正弦波和构成的信号:

L = 32;帧尺寸%Fs = 128;%抽样率res = Fs/L;频率分辨率%F1 = 40;第一个正弦波频率F2 = f1 + res;第二正弦波频率Sn1 = dsp。正弦波(频率=f1, SampleRate=Fs, SamplesPerFrame=L);Sn2 = dsp。正弦波(频率=f2, SampleRate=Fs, SamplesPerFrame=L，振幅=2);X = sn1() + sn2();

计算的FFTx并绘制FFT的大小。注意两个正弦波在FFT的相邻样本中。这意味着你不能区分更接近的频率 $\frac{{\mathrm{F}}_{\mathit{年代}}}{\mathit{l}}$．

X = fft(X);scopeFFT = dsp。ArrayPlot (SampleIncrement = Fs / L,.．.XOffset = - f / 2,.．.包含=“频率(赫兹)”，.．.YLabel =“级”，.．.Title =“双边谱”，.．.YLimits = [-0.1 - 1.1]);scopeFFT (abs (fftshift (X) / L);

变焦FFT

假设您有一个应用程序，您只对奈奎斯特区间的子带感兴趣。变焦FFT背后的思想是通过在较短的信号上计算小FFT来保持与原始信号上全尺寸FFT相同的分辨率。较短的信号来自于抽取原始信号。节省来自能够计算更短的FFT，同时实现相同的分辨率。这是直观的:对于抽取因子D，新的采样率为 ${\mathit{F}}_{\mathrm{sd}}＝\frac{{\mathrm{F}}_{\mathit{年代}}}{\mathit{D}}$，新的帧大小(和FFT长度)为 ${\mathit{l}}_{\mathit{d}}＝\frac{\mathit{l}}{\mathit{D}}$，则抽取信号的分辨率为 $\frac{{\mathit{F}}_{\mathrm{sd}}}{{\mathit{l}}_{\mathit{d}}}＝\frac{{\mathrm{F}}_{\mathit{年代}}}{\mathit{l}}$．

DSP系统工具箱提供放大FFT功能的MATLAB和Simulink，通过dsp。ZoomFFT系统对象和缩放FFT库块。下一节将更详细地讨论放大FFT算法。

混合器方法

在讨论算法使用之前dsp。ZoomFFT，我们提出了混合方法，这是一种流行的放大FFT方法。

对于这里的示例，假设您只对间隔[16 Hz, 48 Hz]感兴趣。

BWOfInterest = 48 - 16;Fc = (16 + 48)/2;感兴趣带宽的中心频率

可实现的抽取因子为:

BWFactor =楼层(f /BWOfInterest)

BWFactor = 4

混频器方法包括首先使用混频器将感兴趣的频段移到直流，然后进行低通滤波和抽取BWFactor．

设计抽取滤波器的系数designMultirateFIR函数。的dsp。FIRDecimator系统对象采用高效的多相FIR抽取结构实现低通和降采样。

B = designMultirateFIR(1, BWFactor);D = dsp。FIRDecimator(BWFactor, B);

现在，将信号混合到直流，并通过FIR抽取器对其进行滤波:

运行几个输入帧以消除FIR瞬态响应为K = 1:10抓取输入信号的一帧X = sn1()+sn2();混音到DCindVect =(0:数字(x)-1)。' + (k-1) * size(x,1);y = x .* exp(-2*pi*indVect*Fc*1j/Fs);通过FIR抽取器进行过滤xd = D(y);结束

现在取滤波信号的FFT(注意，与常规FFT相比，FFT长度减少了BWFactor，或抽取长度，同时保持相同的分辨率):

Xd = fft(Xd);Ld = L/BWFactor;Fsd = Fs/BWFactor;scopeMixing = dsp。ArrayPlot (SampleIncrement = Fs / L,.．.XOffset = Fsd / 2,.．.包含=“频率(赫兹)”，.．.YLabel =“级”，.．.Title =放大FFT频谱:混频器方法，.．.YLimits = [-0.1 - 1.1]);scopeMixing (abs (fftshift (Xd)) / Ld);

复值混合器为每个高速率样本增加一个额外的乘法，这是不有效的。下一节将介绍另一种更有效的放大FFT方法。

带通采样

另一种放大FFT方法利用了带通滤波(有时也称为欠采样)的已知结果:假设我们对带感兴趣 $［{\mathit{F}}_1，{\mathit{F}}_2］，$一个具有采样率的信号 ${\mathit{F}}_{\mathit{年代}}$．如果我们让信号通过一个复杂(单侧)带通滤波器，中心在 ${\mathit{F}}_{\mathit{c}}＝\frac{{\mathit{F}}_1+{\mathit{F}}_2}{2}$还有带宽 $\mathrm{BW}＝{\mathit{F}}_2-{\mathit{F}}_1$，然后将其采样降低到 $D＝\lfloor \frac{{\mathit{F}}_{\mathit{年代}}}{\mathrm{BW}}\rfloor$，我们将把所需的波段下拉到基带。

一般来说，如果 ${\mathit{F}}_{\mathit{c}}$不能用形式表示 $k\frac{{\mathit{F}}_{\mathit{年代}}}{\mathit{D}}$(其中K为整数)，则移位后抽取的频谱将不会以DC为中心。事实上，中心频率Fc将被转换为[2]:

在这种情况下，我们可以使用混频器(在抽取信号的低采样率下运行)将所需的频段集中到零赫兹。

使用上一节的例子，我们通过调制所设计的低通滤波器的系数来得到复带通滤波器的系数:

N =长度(d .分子);Bbp = B .*exp(1j*(Fc / Fs)*2*pi*(0:N-1));Dbp = dsp。FIRDecimator(BWFactor, Bbp);

现在执行过滤和FFT:

为K = 1:10运行几次以消除过滤器中的瞬态X = sn1()+sn2();xd = Dbp(x);结束Xd = fft(Xd);scopeBandpassSampling = dsp。ArrayPlot (SampleIncrement = Fs / L,.．.XOffset = Fsd / 2,.．.包含=“频率(赫兹)”，.．.YLabel =“级”，.．.Title =放大FFT频谱:带通采样方法，.．.YLimits = [-0.1 - 1.1]);scopeBandpassSampling (abs (fftshift (Xd)) / Ld);

使用多速率，多级带通滤波器

在前一节中使用的FIR抽取器是一个单级多速率滤波器。我们可以通过使用多级设计来减少滤波器的计算复杂性(事实上，这是在dsp。ZoomFFT)．看到多级汇率转换欲知详情。

考虑下面的例子，其中输入采样率是48 KHz，感兴趣的带宽是区间[1500,2500]Hz。可实现的抽取因子是 $\lfloor \frac{48000}{1000}\rfloor ＝48$．

首先，设计一个单级decimator:

Fs = 48e3;Fc = 2000;%带通滤波器中心频率BW = 1e3;感兴趣的带宽Ast = 80;阻带衰减%P = 12;多相长度%BWFactor =楼层(Fs/BW);B = designMultirateFIR(1,BWFactor,P,Ast);N =长度(B);Bbp = B .*exp(1j*(Fc / Fs)*2*pi*(0:N-1));D_single_stage = dsp。FIRDecimator（BWFactor，Bbp);

现在，使用多级设计实现相同的抽取器，同时保持与单级情况相同的阻带衰减和转换带宽(参见kaiserord关于过渡宽度计算的详细信息):

tw = (Ast - 7.95) / (N * 2.285);D_multi_stage = designMultistageDecimator(BWFactor, Fs, tw*Fs/(2*pi)， Ast);流(“过滤级数:%d\n”， getNumStages(D_multi_stage));

过滤级数:5级

为ns = 1: D_multi_stage。getNumStages stgn = D_multi_stage。（“舞台”+ ns);流(“阶段%i:抽取因子= %d, FIR长度= %d\n”，.．.ns, stgn。DecimationFactor,.．.长度(stgn.Numerator));结束

第一阶段:抽取因子= 2,FIR长度= 7第二阶段:抽取因子= 2,FIR长度= 7第三阶段:抽取因子= 2,FIR长度= 11第四阶段:抽取因子= 2,FIR长度= 15第五阶段:抽取因子= 3,FIR长度= 75

注意，D_multi_stage是一个五级多速率低通滤波器。我们将其转换为带通滤波器，通过对每个级的系数进行频移，同时考虑到累积抽取因子:

Mn = 1;累积抽取因子进入阶段n为ns = 1: D_multi_stage。getNumStages stgn = D_multi_stage。（“舞台”+ ns);num = stgn.分子;N =长度(num);num = num .*exp(1j*(Fc * Mn/ Fs)*2*pi*(0:N-1));stgn。分子= num;Mn = Mn*stgn.DecimationFactor;结束

比较单级滤波器和多级滤波器的成本，后者的计算效率明显更高。

对于单级过滤器，此成本为:

成本(D_single_stage)

ans =带字段的结构:NumCoefficients: 1129 NumStates: 1104 MultiplicationsPerInputSample: 23.5208 AdditionsPerInputSample: 23.5000

而多级滤波器的成本为:

成本(D_multi_stage)

ans =带字段的结构:NumCoefficients: 77 NumStates: 108 MultiplicationsPerInputSample: 6.2500 AdditionsPerInputSample: 5.2917

接下来，比较两个滤波器的频率响应，并验证它们具有相同的通带特性。阻带的差异可以忽略不计。

vi = fvtool(D_single_stage,D_multi_stage,DesignMask=“关闭”传说=“上”）;传奇(vis、“单级”，“多级”）

最后利用多级滤波器进行放大FFT:

Fftlen = 32;L = BWFactor * fftlen;音调= [1625 2000 2125];%正弦波音Sn = dsp。正弦波(SampleRate=Fs, Frequency=tones, SamplesPerFrame=L);Fsd = Fs / BWFactor;%计算FFT的频率点F = Fc + Fsd/fftlen*(0:fftlen-1)-Fsd/2;步进带通decimator为k=1:100 x = sum(sn()，2) + 1e-2 * randn(L,1);y = D_multi_stage(x);结束绘制光谱输出scopeZoomFFT = dsp。ArrayPlot (XDataMode =“自定义”，.．.CustomXData = F,.．.YLabel =“级”，.．.包含=“频率(赫兹)”，.．.YLimits = [0, 1],.．.标题= sprintf (“变焦FFT:分辨率= %f Hz”, (Fs / BWFactor) / fftlen));Z = fft(y,fftlen,1);Z = fftshift(Z);(abs(z)/fftlen) release(scopeZoomFFT)

使用dsp。ZoomFFT

dsp。ZoomFFT是一个系统对象，它实现了基于前一节中描述的多速率多级带通滤波器的放大FFT。您指定所需的中心频率和抽取因子，以及dsp。ZoomFFT将设计滤波器并将其应用于输入信号。

使用dsp。ZoomFFT放大前一节示例中的正弦音调:

zfft = dsp.ZoomFFT(BWFactor,Fc,Fs);为k=1:100 x = sum(sn()，2) + 1e-2 * randn(L,1);Z = zfft(x);结束Z = fftshift(Z);(abs(z)/fftlen) release(scopeZoomFFT)