Документ взят из кэша поисковой машины. Адрес оригинального документа : http://lnfm1.sai.msu.ru/grav/russian/lecture/filtr/labs/Lab3arss.m Дата изменения: Thu Feb 24 17:52:11 2011 Дата индексирования: Mon Oct 1 23:42:10 2012 Кодировка:

% program is written 20.01.2009 by L.V. Zotov
clear;

N_signal=1024;
% generating two-sin signal
for (k=1:1:N_signal)
garm(k)=0.05*sin(2*pi/10*(k-1))+0.05*sin(2*pi/100*(k-1));
end;

plot(garm);

% ARMA process generating
ar(1)=0.1*randn(1);
ar(2)=-0.2*ar(1)+0.1*randn(1);
for (i=3:1:N_signal)
ar(i)=0.7*ar(i-1)-0.2*ar(i-2)+0.1*randn(1);
end;

plot(ar);

%making a sum
signal=ar;%garm+

plot(signal);
%--------------------------------------------------------------------------
%-------------

%output to the file
%chnage the path
cd /home/leonid/lectures/labs/lab2;
fout=fopen('signal.dat','wt');
for(i=1:1:N_signal)
fprintf(fout,'%12.6f%12.6f\n',i,signal(i));
end;
fclose(fout);

savefile = 'signal.mat';
save(savefile,'signal');

% fast Fourier transformation
N_ft=N_signal;

Ftrns_signal=fft(signal,N_ft);

% frequency calculation
% for the half of the spectrum not counting 1 coefficient - an avarage
% N_ft is odd or even - ?
dt=1;
N_half=round((N_ft-1)/2);
freq=(1:N_half)/N_ft/dt;

ampl_spectrum=abs(Ftrns_signal)/N_ft;

energy_spectrum=Ftrns_signal.*conj(Ftrns_signal)/N_ft^2;

% plot half of the spectrum without multiplication by 2
subplot(2,1,1);
plot(freq, ampl_spectrum(2:N_half+1));
title('amplitude spectrum ')
xlabel('frequency')

subplot(2,1,2);
plot(freq, energy_spectrum(2:N_half+1));
title('energy spectrum ')
xlabel('frequency')

clf;

%--------------------------------------------------------------------------------------------------
% avarage substraction
abs(Ftrns_signal(1))/N_ft %equal
mean(signal)
signal_centered=signal-mean(signal);

% ACF calculation

for(tau=1:1:N_signal)
acf(tau)=0;
for(j=1:1:N_signal-tau)
acf(tau)=acf(tau)+signal_centered(j)*signal_centered(j+tau-1);
end;
acf(tau)= acf(tau)/(N_signal-tau+1);
end;
plot(acf);

plot(acf(1:254));

%--------------------------------------------------------------------------------------------------

% Power Spectrum Dencity calculations

Ftrns_acf=fft(acf,N_ft);

plot(abs(Ftrns_acf)/N_ft);

power_spectrum=abs(Ftrns_acf)/N_ft;
plot(freq, power_spectrum(2:N_half+1));

%--------------------------------------------------------------------------------------------------


% Window generating
% half of simmetric window is needed for positive frequencies
N_half_window=254
for(i=1:1:N_half_window)
half_window(i)=(1-abs(i-1)/(N_half_window-1));
end;

plot(-N_half_window:-1,half_window(N_half_window:-1:1),1:N_half_window,half_window)

% ACF by window multiplication
for(i=1:1:N_half_window)
acf_w(i)=acf(i)*half_window(i);
end;

plot(acf_w);

%--------------------------------------------------------------------------------------------------

% Blackman-Tukey power spectrum calculation

Ftrns_acf_w=fft(acf_w,N_ft);

power_spectrum_w=abs(Ftrns_acf_w)/N_ft;
plot(freq, power_spectrum_w(2:N_half+1));