このページの翻訳は最新ではありません。ここをクリックして、英語の最新版を参照してください。
ズーム 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)
ズーム FFT Simulink ブロック
ズーム FFT ブロックは dsp.ZoomFFT の機能を Simulink で利用できるようにします。モデル dspzoomfft で、44100 Hz でサンプリングされた入力信号の周波数帯 [800 Hz, 1600 Hz] を、ズーム FFT ブロックを使用して検査します。
参考文献
[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)