Wolfram 信号処理 | 試してみよう

コードはどれも，コード内をクリックしてから，編集したり，
+
を押して実行したりできます．
信号の処理と解析 包括的な信号処理機能は，強力な微積分，統計，機械学習の機能に強固に統合され，工学，金融学，医学等の学生およびプロフェッショナルに最適なものとなっています．

ローパスIIRバターワースフィルタ

サンプリング周期が
T=1
と仮定して，カットオフ周波数が
ω
c
=π/6
の無限インパルス応答(IIR)フィルタで次数が
n=3
のアナログバターワースプロトタイプフィルタを作成する：
実行
In[]:=
T=1;
ω
c
=
π
6
;n=3;​​
Ω
c
=
2.
T
Tan
ω
c
2
;​​tf=ButterworthFilterModel[{n,
Ω
c
},s]//TransferFunctionExpand//Chop
双一次変換を使って，対応するデジタルフィルタを作成する：
実行
In[]:=
dtf=ToDiscreteTimeModel[tf,T,z]//Chop
このフィルタのボード線図を描く：
実行
In[]:=
BodePlotdtf,{0,π},


IIRフィルタを使ったデータのフィルタリング

ノイズの多い正弦波信号を作成する：
実行
In[]:=
sig=TableSin
3.π
128
n,{n,0,127}+RandomReal[{-0.25,0.25},128];
IIRバターワースフィルタを定義する：
実行
In[]:=
dtf=
0.153903
3
(1.+z)
-4.70766+19.0258z-26.6767
2
z
+13.5898
3
z
1
;
ノイズの多い信号にフィルタを掛ける：
実行
In[]:=
res=RecurrenceFilter[dtf,sig,Padding->None];​​ListPlot[{sig,res}]

窓関数法を使ったローパスFIRフィルタ

カットオフ周波数が
ω
c
=π/6
で長さが
L=17
の有限インパルス応答(FIR)フィルタを作成する：
実行
In[]:=
ω
c
=
π
6
;L=17;​​h=LeastSquaresFilterKernel[{"Lowpass",
ω
c
},L]
単位サンプル応答にハン窓を適用する：
実行
In[]:=
w=Array[HannWindow,L,{-0.5,0.5}];​​fir=wh
結果のフィルタの周波数応答を表示する：
実行
In[]:=
PlotEvaluate[20Log10[Abs[ListFourierSequenceTransform[fir,ω]]]],{ω,0,π},


FIRを使ったデータのフィルタリング

ノイズの多い信号を作成する：
実行
In[]:=
sig=Table[SquareWave[n/64],{n,0,127}]+RandomReal[{-0.25,0.25},128];
ローパスFIRフィルタを作成する：
実行
In[]:=
fir=LeastSquaresFilterKernel[{"Lowpass",π/6},17]Array[HannWindow,17,{-0.5,0.5}];
ListConvolve
を使って信号にフィルタを掛ける：
実行
In[]:=
res=ListConvolve[fir,sig,9];​​ListPlot[{sig,res}]
同様に，
LowpassFilter
関数を使って信号にフィルタを掛ける：
実行
In[]:=
ListPlot[{sig,LowpassFilter[sig,π/6,17,HannWindow,Padding->0]}]

DTMF信号のパワースペクトル

ダイアルトーンの数字2について，デュアルトーン多重周波数(DTMF)で正規化された一対の周波数（単位：ラジアン/サンプル）を定義する：
実行
In[]:=
{fp,fs}={697,1336};​​{wp,ws}=
2π
8000.
fp,
2π
8000.
fs
約200ミリ秒の間のデュアルトーン信号を作成する：
実行
In[]:=
dtmf=Table[Sin[wpn]+Sin[wsn],{n,0,1799}]+RandomReal[{-0.25,0.25},1800];​​audio=AudioPad[Audio[dtmf,SampleRate8000],{0.02,0.02}]
パワースペクトルを表示する：
実行
In[]:=
Periodogramaudio,

ウェルチの方法は，平滑化された重なり合うパーティションのパワースペクトルを平均化する：
実行
In[]:=
Periodogramaudio,600,300,BlackmanWindow,


スペクトログラムを使った信号解析

ダイアルトーンの数字9と1について，DTMFで正規化された対の周波数（単位：ラジアン/サンプル）を定義する：
実行
In[]:=
freqs=
2π
8000.
#[[1]],
2π
8000.
#[[2]]&/@{{852,1447},{697,1209}}
デュアルトーンシーケンスの音声クリップを作成する：
実行
In[]:=
dtmf=(Table[Sin[#[[1]]n]+Sin[#[[2]]n],{n,0,1799}]+RandomReal[{-0.25,0.25},1800]&)/@freqs;​​audio=AudioJoin[AudioPad[Audio[#,SampleRate8000],{0.01,0.01}]&/@dtmf]
スペクトログラムを見る：
実行
In[]:=
Spectrogramaudio,

セグメントの長さを長くすることで周波数分解能を向上させる：
実行
In[]:=
Spectrogramaudio,512,

さらに，セグメント間のオーバーラップを大きくすることで時間分解能も向上させる：
実行
In[]:=
Spectrogramaudio,512,64,


時系列処理

日々の気温データを得る．データのサンプリング期間は
T=1
日（つまり
T=86400
秒）であり，サンプルレートは1秒間につき
sr=
1
86400
となる：
実行
In[]:=
ts=TimeSeries
Time:
01 Jan 2024
to
06 Aug 2024
Data points: 219
;
カットオフ周波数を1週間の平均気温に近付けたローパスフィルタを適用することで，データを平滑化する：
実行
In[]:=
sr=1/86400;​​DateListPlotts,LowpassFilterts,
π
6
sr,15,

さらに，おおよその月平均気温を計算して可視化する：
実行
In[]:=
DateListPlotts,LowpassFilterts,
π
6
sr,15,LowpassFilterts,
π
186
sr,65,
