(50) スペクトル密度

波形とスペクトル密度の関係を見ます。通常は、パワースペクトルを見ますが、振幅の小さなところの特徴も見るため、振幅スペクトルを見ます。

まず、振幅スペクトルを計算する関数とそれをプロットする関数を用意しましょう。

/* 8.2.9m */
load("fft")$

matchdeclare(i, integerp)$
defrule(exponent_of_2, 2^i)$

ampSpectrum(x) := block([len, y],
  if exponent_of_2(factor(length(x))) = false then error(length(x), "is not power of 2"),
  len : 1 + length(x)/2,
  y : rest(inverse_fft(x), - (length(x)/2 - 1)),
  y : abs(y),
  y : y/lmax(y),
  return(y)
  )$

spectrumPlot(x) := block([s],
  s : ampSpectrum(x),
  plot2d([discrete, makelist(i, i, 1, length(s)), s])
  )$

使用した機能

rest
lmax

高速Fourier変換(FFT)を利用するために、データ長は2のべき乗以外はエラーとしています。振幅スペクトル密度は、見やすさのため、最大振幅を1に規格化しています。

では、矩形波から見てましょう。

/* 8.2.10m */
sig(j) := 2*mod(floor(j/32), 2) - 1$
x : makelist(sig(j), j, 1, 256);
plot2d([discrete, makelist(i, i, 1, length(x)), x]);
spectrumPlot(x);

元の波形が奇函数なので、奇数の高調波にスパイクが現れます。

二つ目の例は、トーンが途中で消えるような波形です。

/* 8.2.11m */
sig(j) := if j < 45 then sin(2*%pi*j/5) else 0$
x : makelist(sig(j), j, 1, 256)$
plot2d([discrete, makelist(i, i, 1, length(x)), x])$
spectrumPlot(x)$

トーンの周波数にピークが見えますが、全体に小さな成分が波打っています。これは、Fourier変換が波形全体を周期函数と見る結果、x=最大とx=0の間でも意図せず波形に不連続性を持ち込んでいるためです。

このずれを消すためには、Fourier変換の前に窓函数を適用します。
窓函数にHann函数 $\sin^2x$ を使ってみましょう。パワースペクトルではなく、振幅スペクトルを見ているので、 $|\sin x|$ を掛けることになります。

/* 8.2.12m */
hannW(j, n) := sin(%pi*(j - 1)/(n - 1))^2$
x : makelist(float(sig(j) * hannW(j, 45)), j, 1, 256)$
plot2d([discrete, makelist(i, i, 1, length(x)), x])$
spectrumPlot(x)$

スペクトル密度は、一見わかりにくい信号から周波数を取り出すのに有効です。
以下のAM変調された信号を例に見てみましょう。

$x(t)=(1+\frac{1}{10}\cos\frac{2\pi t}{30})\sin\frac{2\pi t}{15}$

/* 8.2.14m */
sig(j) := sin(2*%pi*j/15) * (1 + 0.1*cos(2*%pi*j/30))$
x : makelist(sig(j), j, 1, 512)$
plot2d([discrete, makelist(i, i, 1, length(x)), x], [style, dots])$
x : ampSpectrum(x)$
x : 20*log(x)/log(10)$
plot2d([discrete, makelist(i, i, 1, length(x)), x], [y, -60, 0])$