缩放 FFT

打开实时脚本

此示例说明缩放 FFT，这是一种用于在高分辨率下分析部分频谱的信号处理方法。DSP System Toolbox™ 在 MATLAB® 中通过 dsp.ZoomFFT System object™ 提供此功能，在 Simulink® 中通过 DSP System Toolbox 库中的 Zoom FFT 模块提供此功能。

标准 DFT/FFT 的限制

数字信号的长度确定其频谱分辨率。以下示例对此进行了说明。首先创建一个由两个正弦波之和组成的信号。

L   = 32;       % Frame size
Fs  = 128;      % Sample rate
res = Fs/L;     % Frequency resolution
f1  = 40;       % First sine wave frequency
f2  = f1 + res; % Second sine wave frequency
 
sn1 = dsp.SineWave(Frequency=f1, SampleRate=Fs, SamplesPerFrame=L);
sn2 = dsp.SineWave(Frequency=f2, SampleRate=Fs, SamplesPerFrame=L, Amplitude=2);

x = sn1() + sn2();

计算 x 的 FFT 并绘制 FFT 的幅值。两个正弦波位于 FFT 的相邻采样中，这意味着您无法区分比 $\frac{F_{s}}{L}$ 更接近的频率。

X = fft(x);
scopeFFT = dsp.ArrayPlot(SampleIncrement=Fs/L, ...
    XOffset=-Fs/2, ...
    XLabel="Frequency (Hz)", ...
    YLabel="Magnitude", ...
    Title="Two-sided spectrum", ...
    YLimits=[-0.1 1.1]);
scopeFFT(abs(fftshift(X))/L);

缩放 FFT

假设您有一项应用，在其中您只对奈奎斯特区间的某个子带感兴趣。缩放 FFT 旨在通过在较短信号上计算小 FFT 来保持您对原始信号使用全大小 FFT 所能达到的相同分辨率。较短的信号来自对原始信号进行抽取的结果。好处在于能够计算更短的 FFT，同时达到相同的分辨率。对于抽取因子 D，新采样率为 $F_{sd} = \frac{F_{s}}{D}$ ，新帧大小（和 FFT 长度）为 $L_{d} = \frac{L}{D}$ ，因此所抽取信号的分辨率为 $\frac{F_{sd}}{L_{d}} = \frac{F_{s}}{L}$ 。

混频器方法

混频器方法是一种常见的缩放 FFT 方法。

假设您只对 [16 Hz, 48 Hz] 的带宽区间感兴趣，并计算该区间的中心频率。

BWOfInterest = 48 - 16;
Fc           = (16 + 48)/2; % Center frequency of bandwidth of interest

使用前一节示例中的采样率计算可实现的抽取因子。

BWFactor = floor(Fs/BWOfInterest)

BWFactor = 
4

混频器方法包括先使用混频器将感兴趣的频带下移到 DC，然后进行低通滤波并按因子 BWFactor 抽取。

使用 designMultirateFIR 函数设计抽取滤波器的系数。返回的 dsp.FIRDecimator System object 使用高效的多相 FIR 抽取结构实现低通和下采样操作。

D = designMultirateFIR(DecimationFactor=BWFactor,SystemObject=true);

将信号下混频到 DC，并通过 FIR 抽取器对其进行滤波。

% Run several input frames to eliminate the FIR transient response
for k = 1:10
    % Grab a frame of the input signal
    x = sn1()+sn2();

    % Downmix to DC
    indVect = (0:numel(x)-1).' + (k-1) * size(x,1);
    y = x .* exp(-2*pi*indVect*Fc*1j/Fs);

    % Filter through FIR decimator
    xd = D(y);    
end

计算滤波后信号的 FFT。为了保持相同的分辨率，与常规 FFT 变换中的 FFT 长度相比，混频器方法中的 FFT 长度会按 BWFactor（即抽取长度）减少。

Xd  = fft(xd);
Ld  = L/BWFactor;
Fsd = Fs/BWFactor;

scopeMixing = dsp.ArrayPlot(SampleIncrement=Fs/L, ...
    XOffset=Fsd/2, ...
    XLabel="Frequency (Hz)", ...
    YLabel="Magnitude", ...
    Title="Zoom FFT Spectrum: Mixer Approach", ...
    YLimits=[-0.1 1.1]);

scopeMixing(abs(fftshift(Xd))/Ld);

复数值混频器为每个高速率采样增加一次额外的乘法运算，这种方式效率不高。下一节将介绍一种更高效的替代缩放 FFT 方法。

带通采样

另一种缩放 FFT 方法利用带通滤波（有时也称为欠采样）的已知结果：假设您对采样率为 $F_{s}$ 的信号的频带 $[F_{1}, F_{2}]$ 感兴趣。如果您使信号通过以 $F_{c} = \frac{F_{1} + F_{2}}{2}$ 为中心、带宽为 $BW = F_{2} - F_{1}$ 的复数（单边）带通滤波器，然后按因子 $D = ⌊ \frac{F_{s}}{BW} ⌋$ 对其进行下采样，则会将所需的频带下移到基带。

通常，如果您无法将 $F_{c}$ 表示为 $k \frac{F_{s}}{D}$ 形式（其中 k 是整数），则经过偏移和抽取的频谱将不会以 DC 为中心。实际上，中心频率 Fc 会转换为 [2]：

$F_{d} = F_{c} - \frac{F_{s}}{D} ⌊ \frac{D \cdot F_{c} + F_{s} / 2}{F_{s}} ⌋$

在这种情况下，您可以使用混频器（以所抽取信号的低采样率运行）使所需的频带以零赫兹为中心。

使用前一节的示例，通过调制所设计的低通滤波器的系数来获得复数带通滤波器的系数。

N   = length(D.Numerator);
numbp = D.Numerator .*exp(1j*(Fc / Fs)*2*pi*(0:N-1));
Dbp = dsp.FIRDecimator(BWFactor, numbp);

执行滤波和 FFT 运算。

for k = 1:10
    % Run a few times to eliminate transient in filter
    x = sn1()+sn2();
    xd = Dbp(x);    
end
Xd = fft(xd);
scopeBandpassSampling = dsp.ArrayPlot(SampleIncrement=Fs/L, ...
                            XOffset=Fsd/2, ...
                            XLabel="Frequency (Hz)", ...
                            YLabel="Magnitude", ...
                            Title="Zoom FFT Spectrum: Bandpass Sampling Approach", ...
                            YLimits=[-0.1 1.1]);
scopeBandpassSampling(abs(fftshift(Xd))/Ld);

多速率多级带通滤波器

前一节中的 FIR 抽取器是单级多速率滤波器。您可以通过改用多级设计来降低滤波器的计算复杂度（实际上，这是 dsp.ZoomFFT System object 中采用的方法）。有关详细信息，请参阅Multistage Rate Conversion。

请考虑以下示例，其中输入采样率为 48 KHz，感兴趣的带宽是 [1500,2500] Hz 区间。则可实现的抽取因子为 $⌊ \frac{48000}{1000} ⌋ = 48$ 。

首先，设计一个过渡带宽度为奈奎斯特频率 1% 的单级抽取器。

Fs       = 48e3;
Fc       = 2000; % Bandpass filter center frequency
BW       = 1e3;  % Bandwidth of interest
Ast      = 80;   % Stopband attenuation
BWFactor = floor(Fs/BW);
TW       = 0.01;
B        = designMultirateFIR(DecimationFactor=BWFactor,...
                            TransitionWidth=TW,...
                            StopbandAttenuation=Ast);

N        = length(B);
Bbp      = B.*exp(1j*(Fc / Fs)*2*pi*(0:N-1));
D_single_stage = dsp.FIRDecimator(BWFactor, Bbp);

现在，使用多级设计实现相同的抽取器，同时保持与单级情况相同的阻带衰减和带宽。您可以通过从抽取器滤波器的截止频率中减去过渡带宽度的一半来计算输入带宽。

Fcd = 1/BWFactor;        % Decimation passband cutoff frequency
BWNormalized = Fcd-TW/2; % Normalized input bandwidth
D_multi_stage = designRateConverter(DecimationFactor=BWFactor, ...
                                    InputSampleRate=Fs, ...
                                    Bandwidth = BWNormalized*(Fs/2), ...
                                    StopbandAttenuation=Ast);
fprintf("Number of filter stages: %d\n", getNumStages(D_multi_stage) );

Number of filter stages: 5

for ns=1:D_multi_stage.getNumStages
    stgn = D_multi_stage.("Stage" + ns);
    fprintf("Stage %i: Decimation factor = %d , FIR length = %d\n",...
            ns, stgn.DecimationFactor,...
            length(stgn.Numerator));
end

Stage 1: Decimation factor = 2 , FIR length = 7
Stage 2: Decimation factor = 2 , FIR length = 7
Stage 3: Decimation factor = 2 , FIR length = 11
Stage 4: Decimation factor = 2 , FIR length = 15
Stage 5: Decimation factor = 3 , FIR length = 59

此滤波器是一个五级多速率低通滤波器。通过对每级的系数执行频移，同时考虑累积抽取因子，将其变换为带通滤波器。

Mn = 1; % Cumulative decimation factor entring stage n
for ns=1:D_multi_stage.getNumStages
    stgn = D_multi_stage.("Stage" + ns);
    num  = stgn.Numerator;
    N    = length(num);
    num  = num .*exp(1j*(Fc * Mn/ Fs)*2*pi*(0:N-1));
    stgn.Numerator = num;
    Mn   = Mn*stgn.DecimationFactor;
end

计算单级滤波器的成本。

cost(D_single_stage)

ans = struct with fields:
                  NumCoefficients: 985
                        NumStates: 960
    MultiplicationsPerInputSample: 20.5208
          AdditionsPerInputSample: 20.5000

计算多级滤波器的成本。此滤波器在计算上明显更高效。

cost(D_multi_stage)

ans = struct with fields:
                  NumCoefficients: 67
                        NumStates: 93
    MultiplicationsPerInputSample: 6.0417
          AdditionsPerInputSample: 5.0833

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

FA = filterAnalyzer(D_single_stage,D_multi_stage,FrequencyRange="centered");

最后，使用多级滤波器执行缩放 FFT：

fftlen = 32;
L      = BWFactor * fftlen;
tones  = [1625 2000 2125]; % sine wave tones
sn     = dsp.SineWave(SampleRate=Fs, Frequency=tones, SamplesPerFrame=L);
Fsd    = Fs / BWFactor;
F      = Fc + Fsd/fftlen*(0:fftlen-1)-Fsd/2; % Discrete frequency points

% Step through the bandpass-decimator 
for k=1:100
    x = sum(sn(),2) +  1e-2 * randn(L,1);
    y = D_multi_stage(x);
end 

% Plot the spectral output
scopeZoomFFT = dsp.ArrayPlot(XDataMode="Custom",...
                   CustomXData=F,...
                   YLabel="Magnitude",...
                   XLabel="Frequency (Hz)",...
                   YLimits=[0 1],...
                   Title=sprintf ("Zoom FFT: Resolution = %f Hz",(Fs/BWFactor)/fftlen));

z = fft(y,fftlen,1);
z = fftshift(z);

scopeZoomFFT( abs(z)/fftlen )
release(scopeZoomFFT)

使用 dsp.ZoomFFT

dsp.ZoomFFT System object 基于前一节中所述的多速率多级带通滤波器实现缩放 FFT。当您指定所需的中心频率和抽取因子后，dsp.ZoomFFT 会设计滤波器并将其应用于输入信号。

使用 dsp.ZoomFFT 放大来自前一节中信号的正弦音调。

zfft = dsp.ZoomFFT(BWFactor,Fc,Fs);
for k=1:100
    x = sum(sn(),2) + 1e-2 * randn(L,1);
    z = zfft(x);
end 

z = fftshift(z);
scopeZoomFFT( abs(z)/fftlen)

release(scopeZoomFFT)

参考

[1] Multirate Signal Processing - Harris (Prentice Hall).

[2] Computing Translated Frequencies in digitizing and Downsampling Analog Bandpass - Lyons (https://www.dsprelated.com/showarticle/523.php)

另请参阅

dsp.ZoomFFT | Zoom FFT | designMultirateFIR

缩放 FFT

标准 DFT/FFT 的限制

缩放 FFT

混频器方法

带通采样

多速率多级带通滤波器

使用 dsp.ZoomFFT

参考

另请参阅

主题