2

「buttord」および「butter」関数を使用してバターワース係数を計算するのが困難です。私の目的は、作成した時系列からノイズを除去することです。時系列には、レッド ノイズ成分と周波数 0.3 Hz の正弦波信号があります。時系列のサンプリング周波数は 10 Hz です。

「buttord」のドキュメントhttp://www.mathworks.com/help/signal/ref/buttord.htmlに従って、仕様の [n, Wn] を計算しました (ドキュメントの例 1 に従いました)。

Wp = 0.33/(sfreq/2); Ws = 0.37/(sfreq/2);
passripp = 0.1; stopatten = 40;
[n,Wn] = buttord(Wp,Ws,passripp,stopatten);
[b,a] = butter(n,Wn);
y_butter = filter(b,a,timeseries(:,2));

y_butter を時間でプロットすると、どこでもゼロになります!

'freqz' を使用してフィルターの周波数応答を調べようとしました (512 サンプルを使用)。

freqz(b,a,512,sfreq)

このプロットは、遷移帯域が 1 ~ 4 Hz であることを示しています。

フィルターの背後にある私の理解は次のとおりです。

  • 0.3Hzの信号
  • >> 0.3 Hz のノイズ
  • 0 から 0.33 Hz までのすべてを渡す
  • 0.36 Hz 以上のすべてを減衰

あなたの助けは大歓迎です!

データはここからダウンロードできます: http://dl.dropbox.com/u/1918592/detrendedTS.mat「ts」の列 1 は時間変数、列 2 はデータ変数です

線形近似 (Matlab 'detrend') を detrending して、ウォークアウェイ レッド ノイズ動作の一部を削除しました。

4

1 に答える 1

3

投稿したコードは、57 次 (!!!!) のフィルターを生成して、実行したい操作を実行します。このfreqzコマンドは実際にこのフィルターが有効であることを示していますが、最初に を使用して 2 次セクションに変換せずにそのような高次フィルターを直接実装しようとすると、tf2sos間違いなく重大な数値エラーが発生します。これがおそらく、ゼロしか表示されない理由です。フィルタ係数を 2 次セクションに変換してからフィルタをカスケード接続すると、freqzコマンドを使用して観察されたフィルタ出力が実際に得られるはずです。

そもそもフィルタ次数がこれほど高いのは、開始帯域と停止帯域を非常に近くに配置したためです。buttord単純なローパス バターワース フィルターの場合、この関数を使用する必要はまったくありません。ただ使う...

[b,a] = butter(n,cut/(sfreq/2));

...n希望するロールオフを選択するだけです(6、12、18 dBなど)。倍精度データ表現を使用する場合、約 10 になると問題が発生するだけnですが、比較的小さなフィルター次数で目的の作業が行われることがわかるでしょう。

于 2012-11-08T10:42:12.787 に答える