今日の例題（その１）　 「 ノンパラメトリックなスペクトル推定 」

• 解説：　直接法 と Blackman-Tukey 法　

  演習に先立って解説したように，ノンパラメトリックなスペクトル推定の方法には，信号のフーリエ変換をもとに次の定義式
  

  にしたがって直接計算する方法（ただし，極限を取る操作も期待値を取る操作も無視）と，自己共分散関数のフーリエ変換がパワースペクトルであるという Wiener-Khintchine の定理
  を使う Blackman-Tukey 法 がある。なお，フーリエ変換には FFT が用いられる。
  

• 例題プログラム：　[ spectrum.m ]

  次のプログラムは，信号 x のパワースペクトルを，直接法とB-T法のそれぞれを使って求めるプログラムである。
  コメントつきのプログラムは こちら

  %
% スペクトル推定（直接法とBlackman-Tukey法）
%
Fs = 100;
T = 2;
dt = 1/Fs;
N = T/dt;
t = (0:dt:T-dt);
x= sin(2*pi*15*t);
subplot(3,2,1); plot(t,x)　　
legend('given signal')　　
%
% スペクトルの計算（直接法）
%
y = fft(x)./N;
df = 1/T;
f = (0:df:Fs-df);
Pxxd = abs(y).^2*T;
subplot(3,2,3); stem(f,Pxxd)
legend('power spectrum')
z = ifft(Pxxd).*N;
Rxxd = real(z)./T;
subplot(3,2,5); plot(t,Rxxd)
legend('auto correlation')
sum(x.^2)*dt/T
sum(Pxxd)*df
%
% スペクトルの計算（Blackman-Tukey法）
%
Rxxb = xcorr(x,x,'biased');
tau = (-T+dt:dt:T-dt);
subplot(3,2,4); plot(tau,Rxxb)
legend('auto correlation')
f = (0:df/2:Fs-df);
w = fft(Rxxb)./(2*N-1);
Pxxb = abs(w).*(2*N-1)*dt;
subplot(3,2,6); stem(f,Pxxb)
legend('power spectrum')
sum(Pxxb)*df/2
%


• 演習：　例題プログラムを実行してみよう

  プログラムを編集して実行してみよう。与える信号をいろいろ工夫して，以下のことについて確認しよう。

  • サンプリング周波数とエリアシング
  • 観測時間と周波数分解能
  • 自己共分散関数の推定における 'biased' は本質的か？

今週の課題（その１）

• ある信号のスペクトルを求めること

  ある信号のスペクトルを求めるために，その信号を，サンプリング周波数100Hzで，1秒間観測したデータと，20Hzで，5秒間観測したデータを用意したとする。これらのデータを解析して，元の信号にどういった周波数成分が含まれているか議論しなさい。なお，元の信号には，50Hz以上の周波数は含まれておらず，また，ランダム・ノイズもそれほど大きくないものとする。
  　データはそれぞれ　 data1.txt　および ， data2.txt　という名で用意してありますので，ダウンロードして使ってください。なおこれらはともに，100個の１次元データ x で，MATLABで作成し，以下のようにしてファイルに書き出したものです。

  fid = fopen('data1.txt','w');
fprintf(fid,'%12.8f\n',x);
fclose(fid)

  スペクトルの計算は，上の例題を参考にして，自分なりにプログラムを作成してください。なお，ファイルからのデータの読み込みなどは，その方法をHELP集などで自分で調べて，対処してください。

今日の例題（その２） 「 線形システムの伝達関数を求める 」

• 解説：　クロススペクトルと伝達関数　

  線形システムの伝達関数（周波数応答関数）は，以下の式のとおり，入力 x のパワースペクトル密度関数 Pxx と入出力のクロススペクトル密度関数 Pyx の比として与えられます。
  

  

• 例題プログラム：　[ tranfer.m ]

  次のプログラムは，線形システムの例として，10点の移動平均処理を行うフィルタを取り上げ，その周波数応答関数を求めるプログラムです。適当な２つの正弦波の合成信号をこのフィルタに入力し，そのときの入出力スペクトルから上の式にもとづいて伝達関数を計算するようになっており，その結果を理論的な応答関数と比較しています。
  オリジナルのプログラムは， こちらから，ダウンロード

  %
% 線形システム（移動平均フィルタ）の伝達関数推定
%
Fs = 1000;
T = 1;
t = (0:1/Fs:T);
x = sin(2*pi*50*t) + sin(2*pi*120*t) ; % <----　入力信号（かえてみよう）
%
%　Welch's method によるパワースペクトルの推定
nfft = 256;
window = hanning(256);
noverlap = 128;
dflag = 'none';
[Pxx,f] = psd(x,nfft,Fs,window,noverlap,dflag); %　パワースペクトル
subplot(2,2,1); plot(f,10*log10(Pxx)) 　　　　　% dB表示
legend('Pxx')
%
h = ones(1,10)/10;　　　　　%　10点移動平均フィルタ
Ht = freqz(h,1,f,Fs);　　　%　freqz は,フィルタの理論周波数応答を求める関数
subplot(2,2,2); plot(f,abs(Ht));
　　　%　|A(f)| legend('theoretical response')
%
y = filter(h,1,x);　　　　% filter は，入力 x に対してフィルタからの出力を求める関数
[Pxy,f] = csd(x,y,nfft,Fs,window,noverlap,dflag); % クロススペクトル
[Pyy,f] = psd(y,nfft,Fs,window,noverlap,dflag);
subplot(2,2,3); plot(f,10*log10(abs(Pyy)))
legend('Pyy')
%
Hc = Pxy./Pxx; %　周波数応答の計算
subplot(2,2,4); plot(f,abs(Hc));
legend('estimated response')


• 演習：　例題プログラムを実行してみよう

  まず実行してみて，このプログラムの内容をよく理解してください。
  さて，なぜ上の式で計算した周波数応答が理論値と一致しないのでしょうか？

  • 与える入力信号をいろいろかえてみてください。
  • 入力していない周波数は，当然その応答を推定できるわけがない！
  • 入力にホワイトノイズを加えると，何故うまくいくの？
  • 入力を矩形波にしてみたら，うまくいく？

今週の課題（その２）

• 未知のシステムの伝達関数を求めること

  あるフィルタが，ファンクションＭ-ファイル（関数副プログラムのようなもの）として用意してあります。， こちらから，ダウンロードしてきて，自分のディレクトリにおき，上のプログラムに習って，このフィルタの伝達関数を求めるプログラムを作成してください。ただしこのファンクションＭ-ファイルは，アルゴリズムを隠蔽するため，テキストファイル(unknown_F.m)ではなく，p-code ファイル(unknown_F.p)に変換してあります。
  ファンクションＭ-ファイルをプログラム（スクリプトＭ-ファイル）中で使うのは至って簡単で，一般のMATLAB関数と同じように，そのファンクションＭ-ファイルの名前を定められた引数をつけて引用するだけです。今の場合，入力信号 x を引数とし，そのフィルタ出力を計算する関数として定義していますので，プログラム中では以下のように記述するだけです。

  ｙ = unknown_F(x);

• 課題をレポートとして来週の時間に提出せよ

  ２つの課題それぞれについて，作成したプログラムと結果を印刷し，考察とまとめて，来週の授業時間に提出してください。