clear
t=0:0.001:0.6;                    #サンプリング周波数1000Hz
x=sin(2*pi*100*t)+sin(2*pi*200*t); #100Hzと200Hzの信号
y=x+2*randn(size(t));             #ランダム波と合成
subplot(2,1,1)
grid "on"
title("ランダム波形")
xlabel("(sec)")
plot(t,y)

Y=fft(y,512);                     #フーリエ変換により周波数領域に変換
Pyy=Y.*conj(Y)/512;               #Pyy−スペクトル密度　conj−複素共役
f=1000*(0:255)/512;               #最初の256点を周波数軸上にプロット
subplot(2,1,2)
title("スペクトル")
xlabel("(Hz)")
plot(f,Pyy(1:256))