ローパスフィルターを​用いて、カットオフ周​波数以上の周波数をカ​ットしたい

lowpass関数を用いて、カットオフ周波数以上の周波数をカットしたいのですがうまくいきません。

トンネル電流zに参照信号sin(wt)を掛けたものに、ローパスフィルターを適用して周波数をカットしようとしています。

カットオフ周波数は500Hzに設定していますので、FFTをグラフ化した時に、500Hzより高い周波数のピークが消えるはずなのですが消えません。

pixel_image = 256;     %ラスタ走査によって得られる画像のピクセル数を入力(2^nを入力)
dr = 1/(2*sqrt(3));    %ディザ円半径を入力[格子]
a_fast_grid = 10;      %fast軸走査範囲[格子]
a_slow_grid = 10;      %slow軸走査範囲[格子]
fm=5000;               %ディザ円変調周波数[Hz]
fs= fm*240 ;        %サンプリング周波数[Hz]
f_fast = 10.2;         %走査周波数[Hz]を入力(1[s]の1line走査回数)
start_point_x = 0;     %走査開始点のx座標を入力(1[格子]分動かしたい時は1を入力)
start_point_y = 0;     %走査開始点のy座標を入力(1[格子]分動かしたい時は1を入力)

%fast軸三角波のパラメータ設定
amplitude_fast = a_fast_grid/2;    %fast軸振幅

%slow軸三角波のパラメータ設定
amplitude_slow = a_slow_grid/2;    %slow軸振幅
f_slow = (f_fast)/(2*pixel_image); %slow軸三角波周波数

% 時間ベクトルの生成
total_time=256/f_fast;   %全走査時間
t = linspace(0, total_time, fs * total_time);

x_raster = start_point_x + amplitude_fast*(2/pi)*acos(cos(2*pi*f_fast*t));
y_raster = start_point_y + amplitude_slow*(2/pi)*acos(cos(2*pi*f_slow*t));
x_dither = dr*cos(2*pi*fm*t);
y_dither = dr*sin(2*pi*fm*t);
x = x_raster + x_dither; 
y = y_raster + y_dither;
z1 = cos(2*pi*((x-y)/(sqrt(3))));
z2 = cos(2*pi*(2*y/(sqrt(3))));
z3 = cos(2*pi*((x+y)/(sqrt(3))));
z  = (z1 + z2 + z3);

%参照信号作成
w = 2*pi*fm ;
phi = 0 ;
reference_signal = sin(w*t+phi) ;

%ミキシング信号作成
mixising_signal = z.* reference_signal ;

%ローパスフィルターの適用
f_cutoff = 500 ; %カットオフ周波数[Hz]の設定
lowpass_signal =  lowpass(mixising_signal,f_cutoff,fs) ;


%FFTしてみた結果-------------------------------------------

n = length(t); % サンプル数
f = (0:n-1)*(fs/n); % 周波数軸 [Hz]
half_n = floor(n/2); % データの半分

fft_z = abs(fft(z) / n);
fft_mixising_signal = abs(fft(mixising_signal) / n);
fft_lowpass_signal = abs(fft(lowpass_signal) / n);

figure;
tiledlayout(3,1)
nexttile
plot(f(1:half_n), fft_z(1:half_n)); 
title('トンネル電流のFFT'); xlabel('Frequency [Hz]'); ylabel('Amplitude');
xlim([0 20000]),ylim([0 0.3]);

nexttile
plot(f(1:half_n), fft_mixising_signal(1:half_n)); 
title('ミキシング信号のFFT'); xlabel('Frequency [Hz]'); ylabel('Amplitude');
xlim([0 20000]),ylim([0 0.3]);

nexttile
plot(f(1:half_n), fft_lowpass_signal(1:half_n)); 
title('ローパス適用後のFFT'); xlabel('Frequency [Hz]'); ylabel('Amplitude');
xlim([0 20000]),ylim([0 0.3]);

Answer 1

[Hiroshi Iwamura]

フィルターの遮断周波数／ｆｓが小さいため、フィルターを適切に設計する必要があります。

以下を参照してください。

https://jp.mathworks.com/help/releases/R2024b/signal/ref/lowpass.html#mw_3aa14271-5d7b-4a2a-8b29-b1f47d3ce414

lowpass 関数の Steepness はデフォルトで 0.85なので、

W = (1 - 0.85) * (fs/2 - f_cutoff)

≓ 89925 となってしまいます。

lowpass_signal = lowpass(mixising_signal,f_cutoff,fs,ImpulseResponse="iir",Steepness=0.95) ;

などとするか他の方法でLPFを設計してください。

Steepness=0.95

また、lowpass 関数は戻り値なしで呼び出すとフィルター適用前後の周波数特性を描いてくれるので、それと比べてみると便利かと思います。

lowpass(mixising_signal,f_cutoff,fs,ImpulseResponse="iir",Steepness=0.95) 
subplot(2,1,2)
xlim([0 20])  % デフォルトのX軸単位がkHz

Answer 2

[薫]

ご回答ありがとうございます！

続けて質問させてください。

lowpass_signal = lowpass(mixising_signal,f_cutoff,fs,ImpulseResponse="iir",Steepness=0.95) ;

上記で使っているlowpass関数のローパスフィルター自身の周波数応答を、グラフ化する方法はありますでしょうか？

グラフから以下のことを調べたいので、教えて頂けたら幸いです。

確認したいこと

・カットオフ周波数である500Hzでパワースペクトルが-3dBになっているかどうか。

・どのくらいのdB/decadeで減衰しているか。(すなわち周波数が10倍になったらdBはどのくらい減衰するかを見たい。)

・何次のローパスフィルターなのか。

よろしくお願いいたします。

Answer 3

[Hiroshi Iwamura]

リンク先の公式ドキュメントにあるように

[lowpass_signal d] = lowpass(～

でフィルターオブジェクトが得られ、

>> d

>> filterAnalyzer(d)

等でフィルターの属性や特性が確認できます。

が、最初に書いたようにfsが遮断周波数に対して高すぎるため、フィルター仕様としてはかなり厳しいものとなります。

直接的な設計例は載せますが、フィルター特性も気にするのであれば信号の性質と目的に合わせてちゃんとフィルター設計（マルチレート信号処理、システム全体の見直しを含め）した方が良いかと思います。

デジタル フィルターの設計

マルチレート信号処理

f_cutoff = 500;
fs = 1200000;
fc_normalized = f_cutoff / (fs / 2); % 正規化カットオフ周波数
 
filter_order = 5000;

% 窓関数法 FIR
b = fir1(filter_order, fc_normalized, 'low', hann(filter_order + 1));

fvtool(b, 1, 'Fs', fs);
title('FIR ローパスフィルタの周波数特性');

d = designfilt('lowpassfir', ... 
       'FilterOrder',4095, ... 
       'PassbandFrequency',f_cutoff * 0.8, ...
       'StopbandFrequency',f_cutoff * 1.2, ...
       'DesignMethod','ls', ...
       'PassbandWeight',1, ...
       'StopbandWeight',2, ...
       'SampleRate',fs)
filterAnalyzer(d)

[薫]

ご説明ありがとうございました。

上記の内容を踏まえて、ローパスフィルタを再設計したいと思います。