% spectrum.m 
%
%      スペクトル推定（直接法とBlackman-Tukey法）
%
Fs = 100;                            % サンプリング周波数 (Hz)
T = 1;                               % 観測時間（1sec）　　　　　
dt = 1/Fs;                           % 時間分解能（サンプリング間隔）
N = T/dt;                            % サンプル数
t = (0:dt:T-dt);                     % time vector 
rn = 0.5*randn(size(t));
x = sin(2*pi*2.4*t)+sin(2*pi*30*t)+sin(2*pi*15.4*t)+rn;  % signal      (15Hz)
fid = fopen('data1.txt','w');
fprintf(fid,'%12.8f\n',x);
fclose(fid);
subplot(3,2,1); plot(t,x)
legend('given signal')
%
%      スペクトルの計算（直接法）
%                            
y = fft(x)./N;                       % Discrete Fourier Transform of x
df = 1/T;                            % 周波数分解能
f = (0:df:Fs-df);                    % frequency vector
Pxxd = abs(y).^2*T;                  % パワースペクトル密度 Pxx
subplot(3,2,3); stem(f,Pxxd)
legend('power spectrum')
z = ifft(Pxxd).*N;                   % Inverse Discrete Fourier Transform of Pxx
Rxxd = real(z)./T;                   % パワ－スペクトルのフーリエ逆変換　-＞　自己共分散関数
subplot(3,2,5); plot(t,Rxxd)
legend('auto correlation')
sum(x.^2)*dt/T                       % Avarage Power of the signal
sum(Pxxd)*df
%
%      スペクトルの計算（Blackman-Tukey法）
%
Rxxb = xcorr(x,x,'biased');          % 自己共分散関数 Rxx
tau = (-T+dt:dt:T-dt);               % time(lag) vector
subplot(3,2,4); plot(tau,Rxxb)
legend('auto correlation')
f = (0:df/2:Fs-df);                  % frequency vector
w = fft(Rxxb)./(2*N-1);              % Discrete Fourier Transform of Rxx
Pxxb = abs(w).*(2*N-1)*dt;           % 自己共分散関数のフーリエ変換　-＞　パワースペクトル
subplot(3,2,6); stem(f,Pxxb)
legend('power spectrum')
sum(Pxxb)*df/2                       % Avarage Power of the signal
%
% ---------------------------------------------------------------------------------------
%                                                                2002.5.14  by K.Tsukada