Документ взят из кэша поисковой машины. Адрес оригинального документа : http://lnfm1.sai.msu.ru/grav/english/lecture/filtering/Lab3.m Дата изменения: Tue Jan 27 23:06:56 2009 Дата индексирования: Mon Oct 1 23:38:47 2012 Кодировка: Поисковые слова: photosphere

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

N_signal=1024;
% generating two-sin signal
for (k=1:1:N_signal)
garm(k)=0.5*sin(2*pi/4*(k-1));
end;

plot(garm);

% ARMA process generating
ar(1)=0.5*randn(1);
ar(2)=-0.2*ar(1)+0.5*randn(1);

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

plot(ar);

%making a sum
signal=garm+ar;

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

% fast Fourier transformation
N_1_4=N_signal/4;
N_ft=N_1_4;

Ftrns_signal1=fft(signal(1:N_ft),N_ft);
Ftrns_signal2=fft(signal(N_ft+1:2*N_ft),N_ft);
Ftrns_signal3=fft(signal(2*N_ft+1:3*N_ft),N_ft);
Ftrns_signal4=fft(signal(3*N_ft+1:4*N_ft),N_ft);

Ftrns_signal_avg=(Ftrns_signal1+Ftrns_signal2+Ftrns_signal3+Ftrns_signal4)/4;

% 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;

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

subplot(5,1,2);
plot(freq, abs(Ftrns_signal2(2:N_half+1)));
title('amplitude spectrum 2')
xlabel('frequency')

subplot(5,1,3);
plot(freq, abs(Ftrns_signal3(2:N_half+1)));
title('amplitude spectrum 3')
xlabel('frequency')

subplot(5,1,4);
plot(freq, abs(Ftrns_signal4(2:N_half+1)));
title('amplitude spectrum 4')
xlabel('frequency')

subplot(5,1,5);
plot(freq, abs(Ftrns_signal_avg(2:N_half+1)));
title('amplitude spectrum average')
xlabel('frequency')

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

clf;

% avarage substraction

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));

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


% interferometr

for(i=1:1:N_signal)
if (i<=N_1_4||(i>2*N_1_4&&i<3*N_1_4))
interferometr_window(i)=1;
else
interferometr_window(i)=0;
end;
end;

signal=garm.*interferometr_window;
plot(signal);

N_ft=N_signal;

Ftrns_signal=fft(signal);
Ftrns_itntw=fft(interferometr_window);

% 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;

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

subplot(2,1,2);
plot(freq, abs(Ftrns_itntw(2:N_half+1)));
title('amplitude spectrum interf_window')
xlabel('frequency')



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