Документ взят из кэша поисковой машины. Адрес оригинального документа : http://lnfm1.sai.msu.ru/grav/english/lecture/filtering/Lab2.m Дата изменения: Wed Jan 21 21:17:50 2009 Дата индексирования: Mon Oct 1 23:38:47 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)=sin(2*pi/10*(k-1))+sin(2*pi/100*(k-1));
end;

plot(garm);

% ARMA process generating
ar(1)=0.5*randn(1);

for (i=2:1:N_signal)
ar(i)=-0.9*ar(i-1)+0.5*randn(1);
end;

plot(ar);

%making a sum
signal=garm+ar;

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

%output to the file
%chnage the path
cd C:\work\course\filtr\eng\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(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));