ピークフィルターおよびノッチフィルターの設計

スクリプトを開く

この例では、ピークフィルターおよびノッチフィルターの設計法を示します。特定の周波数でピークに達するまたは刻むフィルターは、信号の特定の周波数コンポーネントを保持または削除する場合に便利です。フィルターの設計パラメーターは、ピークまたはノッチに指定された周波数、および 3-dB の帯域幅またはフィルターの Q ファクターのいずれかです。さらに、指定されたこれらの仕様では、フィルター次数を増加することで、理想フィルターに密接に近似する設計を得ることが可能です。

ノッチフィルター

3000 Hz で抽出された信号の 60 Hz の干渉を排除する必要があるとします。このような目的にはノッチフィルターを使用できます。

F0 = 60;% interference is at 60 HzFs = 3000;% sampling frequency is 3000 Hznotchspec = fdesign.notch('N,F0,Q',2,F0,10,Fs); notchfilt = design(notchspec,'SystemObject',true); fvtool(notchfilt,'Color',“白色”);

フィルターの品質係数または Q ファクターは、希望の周波数における他の周波数からの分離性能を示す尺度です。固定フィルター次数の場合、高次の Q ファクターは極を押してゼロに近づけることで実行されます。

notchspec.Q = 100; notchfilt1 = design(notchspec,'SystemObject',true); fvt= fvtool(notchfilt, notchfilt1,'Color',“白色”); legend(fvt,'Q = 10','Q = 100');

品質係数を指定する他の方法としては、3-dB の帯域幅 BW を指定する方法があります。これらは Q = F0/BW に関連しています。帯域幅を指定する方が、設計されたフィルターの最適な形状を実現するために便利な方法である場合があります。

notchspec = fdesign.notch('N,F0,BW',2,60,5,3000); notchfilt2 = design(notchspec,'SystemObject',true); fvt= fvtool(notchfilt, notchfilt1, notchfilt2,'Color',“白色”); legend(fvt,'Q = 10','Q = 100','BW = 5 Hz');

これまで極を押して安定させることのみが可能だったので、フィルターの急峻な特性をもつ近似を改良するためには、フィルターの次数を増加させる必要があります。

notchspec = fdesign.notch('N,F0,Q',2,.4,100); notchfilt = design(notchspec,'SystemObject',true); notchspec.FilterOrder = 6; notchfilt1 = design(notchspec,'SystemObject',true); fvt= fvtool(notchfilt, notchfilt1,'Color',“白色”); legend(fvt,'2nd-Order Filter','6th-Order Filter');

指定された次数では、通過帯域および/または阻止帯域リップルを可能にすることで鋭い遷移を得ることができます。

N = 8; F0 = 0.4; BW = 0.1; notchspec = fdesign.notch('N,F0,BW',N,F0,BW); notchfilt = design(notchspec,'SystemObject',true); notchspec1 = fdesign.notch('N,F0,BW,Ap,Ast',N,F0,BW,0.5,60); notchfilt1 = design(notchspec1,'SystemObject',true); fvt= fvtool(notchfilt, notchfilt1,'Color',“白色”); legend(fvt,'Maximally Flat 8th-Order Filter',...'8th-Order Filter With Passband/Stopband Ripples',...'Location','NorthWest'); axis([0 1 -90 0.5]);

ピークフィルター

ピークフィルターは、信号から単一の周波数成分 (または周波数の小さな帯域) のみを保持する場合に使用されます。これまでに説明したすべての仕様およびトレードオフは、ピークフィルターに同様に適用されます。例を以下に示します。

N = 6; F0 = 0.7; BW = 0.001; peakspec = fdesign.peak('N,F0,BW',N,F0,BW); peakfilt = design(peakspec,'SystemObject',true); peakspec1 = fdesign.peak('N,F0,BW,Ast',N,F0,BW,80); peakfilt1 = design(peakspec1,'SystemObject',true); fvt= fvtool(peakfilt, peakfilt1,'Color',“白色”); legend(fvt,'Maximally Flat 6th-Order Filter',...'6th-Order Filter With 80 dB Stopband Attenuation','Location','North');

時変ノッチフィルターの実装

時変フィルターを使用するには、シミュレーションの実行中にフィルターの係数を変更しなければなりません。fdesign オブジェクトに基づく自動フィルター設計ワークフローを補完するために、DSP System Toolbox ではその他の機能が提供されており、iirnotch などのフィルター係数を直接計算する関数も含まれています。

動的な (ストリーミングされた) シミュレーション内から呼び出される静的フィルターから始めるとしましょう。この場合、2 次ノッチフィルターが直接作成され、iirnotch で係数が計算されます。設計パラメーターは -3 dB において中心周波数が 1 kHz、帯域幅が 50 Hz、サンプリング周波数は 8 kHz です。

Fs = 8e3;% 8 kHz sampling frequencyF0 = 1e3/(Fs/2);% Notch at 2 kHzBW = 500/(Fs/2);% 500 Hz bandwidth[b, a] = iirnotch(F0, BW)%#okbiquad = dsp.BiquadFilter('SOSMatrix', [b, a]); scope = dsp.SpectrumAnalyzer('PlotAsTwoSidedSpectrum', false,...'SampleRate', Fs,...'SpectralAverages', 16); samplesPerFrame = 256; nFrames = 4096;fork = 1:nFrames x = randn(samplesPerFrame, 1); y = biquad(x); scope(y);end

b = 0.8341-1.1796 0.8341 a = 1.0000 -1.1796 0.6682

設計パラメーター (ノッチフィルターの中心周波数など) は実行時に変化するため、時変フィルターの係数は時間経過とともに変化しなければなりません。設計パラメーターが変化するたびに、係数ベクトル b と a を再計算し、biquad.SOSMatrix を新しい値に設定しなければなりません。これには大量の計算が必要です。このようなタイプのアプリケーションでは、dsp.CoupledAllpassFilter の方が便利なフィルター構造を提供する場合があります。利点としては、本質的な安定性と、係数の設計パラメーターからの分離が挙げられます。

% Build a coupled allpass lattice filter equivalent to biquad[k1, k2] = tf2cl(b, a)%#okapnotch = dsp.CoupledAllpassFilter('Lattice', k1, k2); fvtool(apnotch,'Fs', Fs)

k1 = [] k2 = -0.7071 0.6682

このオールパスベースの構造の利点としては、新しい係数が設計パラメーター F0 および BW から分離されているという点があります。このため、たとえば F0 を 3 kHz にすると、出力は次のようになります。

F0 = 3e3/(Fs/2); [b, a] = iirnotch(F0, BW)%#ok[k1, k2] = tf2cl(b, a)%#ok% NOTE: while a and b both changed, in the lattice allpass form the design% change only affected k2(1)% Now change the bandwidth to 1 kHzBW = 1e3/(Fs/2); [b, a] = iirnotch(F0, BW)%#ok[k1, k2] = tf2cl(b, a)%#ok% The design change now only affected k2(2)% Coefficient decoupling has numerous advantages in real-time systems,% including more economical coefficient update and more predictable% transient behaviour when coefficients change

b = 0.83411.1796 0.8341 a = 1.0000 1.1796 0.6682 k1 = [] k2 = 0.7071 0.6682 b = 0.7071 1.0000 0.7071 a = 1.0000 1.0000 0.4142 k1 = [] k2 = 0.7071 0.4142

以下は上記の原理を適用して、動的なシミュレーション中に設計パラメーターを変更しています。また、フィルター伝達関数の推定に対する影響をライブで可視化しています。実際のアプリケーションでは、一般的に各設計パラメーターを対応するラティスオールパス係数にリンクする個々の式を識別します。主シミュレーションループで iirnotch や tf2cl などのフィルター全体の関数の代わりにこれらの式を使用することで、効率性が向上されます。

% Notch filter parameters - how they vary over timeFs = 8e3;% 8 kHz sampling frequencyf0 = 1 e3 *(0.5、1.5、3、2)/ (Fs / 2);% in [0, 1] rangebw = 500/(Fs/2) * ones(1,4);% in [0, 1] rangemyChangingParams = struct('F0', num2cell(f0),'BW', num2cell(bw)); paramsChangeTimes = [0, 7, 15, 20];% in seconds% Simulation time managementnSamplesPerFrame = 256; tEnd = 30; nSamples = ceil(tEnd * Fs); nFrames = floor(nSamples / nSamplesPerFrame);%对象创建apnotch = dsp.CoupledAllpassFilter('Lattice'); scope = dsp.SpectrumAnalyzer('PlotAsTwoSidedSpectrum', false,...'SampleRate', Fs,...'SpectralAverages', 4,...'RBWSource','Auto'); paramtbl = dspdemo.ParameterTimeTable('Time', paramsChangeTimes,...'Values', myChangingParams,...'SampleRate', Fs/nSamplesPerFrame);% Actual simulation loopforframeIdx = 1:nFrames% Get current F0 and BW[params, update] = paramtbl();if(update)% Recompute filter coefficients if parameters changed[b, a] = iirnotch(params.F0, params.BW); [k1, k2] = tf2cl(b, a);% Set filter coefficients to new valuesapnotch.LatticeCoefficients1{1} = k1; apnotch.LatticeCoefficients2{1} = k2;end% Generate vector of white noise samplesx = randn(nSamplesPerFrame, 1);% Filter noisey = apnotch(x);% Visualize spectrumscope(y);end