ズーム FFT

スクリプトを開く

この例では、ズーム FFT を紹介します。ズーム FFT は、スペクトルの高分解能な解析に使用される信号処理手法です。DSP System Toolbox™ はこの機能を、MATLAB™ では dsp.ZoomFFT System object を介して、Simulink ではズーム FFT ライブラリブロックを介して提供します。

通常の FFT

信号の分解能は長さによって限定されます。これは、簡単な例で示すことができます。2 つの正弦波の和から形成される信号で考えてみましょう。

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 の振幅をプロットしてみます。ここでは 2 つの正弦波が FFT の隣接するサンプルのものであることに注意してください。つまり、Fs/L よりも近い周波数の差を識別することはできません。

X = fft(x);
stem(Fs/L*(0:length(X)-1)-Fs/2,abs(fftshift(X))/L)
grid on;
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('Two-sided spectrum')

ズーム FFT

ナイキスト区間のサブ帯域のみを対象とする応用事例があるとします。ズーム FFT の背景にある概念は、短い信号の小さい FFT を計算して元の信号のフルサイズの FFT から得られる分解能と同じ分解能を得ることです。短い信号は元の信号を間引きしたものです。はるかに短い FFT の計算から同じ分解能が得られれば負荷の軽減になります。直観的には、間引き係数を D とすると、新しいサンプリングレートは Fsd = Fs/D、新しいフレームサイズ (および FFT 長) は Ld = L/D となり、間引き信号の分解能は Fsd/Ld = Fs/L となります。

DSP System Toolbox は MATLAB および Simulink 用のズーム FFT 機能を、dsp.ZoomFFT System object およびズーム FFT ライブラリブロックを介してそれぞれ提供します。この後の節では、ズーム FFT のアルゴリズムについてより詳細に説明します。

混合器手法

dsp.FFT で使用されるアルゴリズムについて説明する前に、一般的なズーム FFT の手法である混合器手法について説明します。

この例では、[16 Hz, 48 Hz] の区間のみを考えるものとします。

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

実現可能な間引き係数を次に示します。

BWFactor = floor(Fs/BWOfInterest);

混合器手法では、最初に混合器を使用して対象の帯域を DC まで下にシフトし、次にローパスフィルター処理と係数 BWFactor による間引きを行います (効率的なポリフェーズ FIR 間引き構造を使用します)。

関数 designMultirateFIR を使用して間引きフィルターの係数を設計します。

B = designMultirateFIR(1,BWFactor);
D = dsp.FIRDecimator('DecimationFactor',BWFactor,'Numerator',B);

ここで信号を混合して DC まで下げ、FIR 間引きによってフィルター処理します。

for k = 1:10
    % Run a few times to eliminate transient in filter
    x = sn1()+sn2();
    % Mix down 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 長は BWFactor つまり間引きの長さによって減少していますが、同じ分解能を維持していることに注意します)。

fftlen = numel(xd);
Xd = fft(xd);
figure
Ld  = L/BWFactor;
Fsd = Fs/BWFactor;
F   = Fc + Fsd/fftlen*(0:fftlen-1)-Fsd/2;
stem(F,abs(fftshift(Xd))/Ld)
grid on
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('Zoom FFT Spectrum. Mixer Approach.')

複素数値の混合器は各高レートサンプルに対して乗算を追加しますが、これは効率的ではありません。次の節では、より効率的な代替手段であるズーム FFT アプローチについて説明します。

バンドパスのサンプリング

代替のズーム FFT 手法はバンドパスフィルター処理の既知の結果を活用します (アンダーサンプリングと呼ばれる場合もあります)。サンプリングレート Fs Hz の信号の帯域 [F1,F2] を考えるものとします。信号を中心が Fc = (F1+F2)/2 で帯域幅が BW = F2 - F1 の複素 (片側) バンドパスフィルターに通し、係数 D = floor(Fs/BW) でダウンサンプリングする場合、目的の帯域はベースバンドまで下がります。

一般的に Fc が k*Fs/D (K は整数) の形式で表現できない場合、シフトした間引きスペクトルの中心は DC になりません。実際には中心周波数 Fc は次のように変換されます [2]。

Fd = Fc - (Fs/D)*floor((D*Fc + Fs/2)/Fs)

この場合、混合器を使用して (間引き信号の低サンプルレートで実行) 目的の帯域の中心を 0 Hz にできます。

前節の例を使用し、設計されたローパスフィルターの係数を変調して複素バンドパスフィルターの係数を取得します。

Bbp     = B .*exp(1j*(Fc / Fs)*2*pi*(0:numel(B)-1));
D = dsp.FIRDecimator('DecimationFactor',BWFactor,'Numerator',Bbp);

フィルター処理と FFT を実施します。

for k = 1:10
    % Run a few times to eliminate transient in filter
    x = sn1()+sn2();
    xd = D(x);
end
Xd = fft(xd);
figure
stem(F,abs(fftshift(Xd))/Ld)
grid on
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('Zoom FFT Spectrum. Bandpass Sampling Approach.')

マルチレート多段バンドパスフィルターの使用

前節のフィルターは単一ステージフィルターです。代わりに多段設計を使用して計算量の少ない複素フィルターを実現できます。これが dsp.ZoomFFT で使用されているアプローチです。

入力サンプルレートが 48 kHz で対象の帯域幅が区間 [1500,2500] Hz である次の例について検討します。実現可能な間引き係数は floor(48000/1000) = 48 です。

最初に単一ステージの間引きを設計します。

Fs       = 48e3;
Fc       = 2000; % Bandpass filter center frequency
BW       = 1e3;  % Bandwidth of interest
Ast      = 80;  % Stopband attenuation
BWFactor = floor(Fs/BW);
B        = designMultirateFIR(1,BWFactor,12,80);
Bbp      = B .*exp(1j*(Fc / Fs)*2*pi*(0:numel(B)-1));
D_single_stage = dsp.FIRDecimator('DecimationFactor',BWFactor,'Numerator',Bbp);

次に多段設計を使用してフィルターを設計しますが、単一ステージの場合と同じ阻止帯域の減衰量と遷移帯域幅を維持します (遷移幅の計算の詳細については、kaiserord を参照してください)。

tw = (Ast - 7.95) / ( numel(B) * 2.285);
fmd = fdesign.decimator(BWFactor,'Nyquist',BWFactor,...
                        'TW,Ast', tw * Fs / (2*pi),Ast,Fs);
D_multi_stage = multistage(fmd,'HalfbandDesignMethod','equiripple','SystemObject',true);
fprintf('Number of filter stages: %d\n', getNumStages(D_multi_stage) );
fprintf('Stage 1: Decimation factor = %d , filter length = %d\n',...
        D_multi_stage.Stage1.DecimationFactor,...
        numel(D_multi_stage.Stage1.Numerator));
fprintf('Stage 2: Decimation factor = %d , filter length = %d\n',...
        D_multi_stage.Stage2.DecimationFactor,...
        numel(D_multi_stage.Stage2.Numerator));
fprintf('Stage 3: Decimation factor = %d , filter length = %d\n',...
        D_multi_stage.Stage3.DecimationFactor,...
        numel(D_multi_stage.Stage3.Numerator));

Number of filter stages: 3
Stage 1: Decimation factor = 6 , filter length = 33
Stage 2: Decimation factor = 2 , filter length = 11
Stage 3: Decimation factor = 4 , filter length = 101

D_multi_stage が 3 ステージのマルチレートローパスフィルターであることに注意してください。これを累積間引き係数を考慮に入れながら各ステージの係数について周波数シフトを実行することによってバンドパスフィルターに変換します。

num1 = D_multi_stage.Stage1.Numerator;
num2 = D_multi_stage.Stage2.Numerator;
num3 = D_multi_stage.Stage3.Numerator;
N1 = length(num1)-1;
N2 = length(num2)-1;
N3 = length(num3)-1;
% Decimation factor between stage 1 & 2:
M12 = D_multi_stage.Stage1.DecimationFactor;
% Decimation factor between stage 2 & 3:
M23 = D_multi_stage.Stage2.DecimationFactor;
num1     = num1 .*exp(1j*(Fc / Fs)*2*pi*(0:N1));
num2     = num2 .*exp(1j*(Fc * M12 / Fs)*2*pi*(0:N2));
num3     = num3 .*exp(1j*(Fc * M12 * M23 / Fs)*2*pi*(0:N3));
D_multi_stage.Stage1.Numerator = num1;
D_multi_stage.Stage2.Numerator = num2;
D_multi_stage.Stage3.Numerator = num3;

単一ステージフィルターと多段フィルターのコストを比較します。

fprintf('Single-stage filter cost:\n\n')
cost(D_single_stage)
fprintf('Multistage filter cost:\n')
cost(D_multi_stage)

Single-stage filter cost:


ans = 

  struct with fields:

                  NumCoefficients: 1129
                        NumStates: 1104
    MultiplicationsPerInputSample: 23.5208
          AdditionsPerInputSample: 23.5000

Multistage filter cost:

ans = 

  struct with fields:

                  NumCoefficients: 113
                        NumStates: 140
    MultiplicationsPerInputSample: 7.0208
          AdditionsPerInputSample: 6.7500

2 つのフィルターの周波数応答も比較します。

vis = fvtool(D_single_stage,D_multi_stage,'DesignMask','off','legend','on');
legend(vis,'Single-stage','Multistage')

多段フィルターを使用してズーム 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;
% Frequency points at which FFT is computed
F      = Fc + Fsd/fftlen*(0:fftlen-1)-Fsd/2;
ap = dsp.ArrayPlot('XDataMode','Custom',...
                   'CustomXData',F,...
                   'YLabel','Magnitude',...
                   'XLabel','Frequency (Hz)',...
                   'YLimits',[0 1],...
                   'Title',sprintf('Zoom FFT. Resolution = %f Hz',(Fs/BWFactor)/fftlen));
for k=1:100
    x = sum(sn(),2) +  1e-2 * randn(L,1);
    y = D_multi_stage(x);
    z = fft(y,fftlen,1);
    z = fftshift(z);
    ap( abs(z)/fftlen )
end
release(ap)

dsp.ZoomFFT の使用

dsp.ZoomFFT は、前節で説明したマルチレート多段バンドパスフィルターに基づいてズーム FFT を実装する System object です。目的の中心周波数と間引き係数を指定すると、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);
    z = fftshift(z);
    ap( abs(z)/fftlen)
end
release(ap)