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

% program is improved 14.02.2011 by L.V. Zotov
clear;


%--------- Part 1---------------------
% fft of artificial signal

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

k=1:1:N_signal;
signal=sin(2*pi/10*(k-1))+sin(2*pi/100*(k-1));

% try to add noise with Coef*randn([1,N_signal])

plot(signal);

%fast Fourier transform calculation

Ftrns_signal=fft(signal);

%amplitude spectrum calculation
apl_spectrum=abs(Ftrns_signal);

plot(apl_spectrum);


%----------- Part 2-----------------------------------------
% download file http://hpiers.obspm.fr/iers/eop/eopc01/eopc01.1900-now
% to the folder (change if needed):
cd C:\work\course\filtr\eng\Lab1;

fin=fopen('eopc01.1900-now','rt');
A=fscanf(fin,'%f',[11 inf]);% A - array of data
fclose(fin);

%determining the size of the signal
l=size(A);
N=l(2);

%selecting the rows of the Array
Date=A(1,1:N);
X_pole=A(2,1:N);
Y_pole=A(4,1:N);
dt=Date(2)-Date(1);

plot3(X_pole,Y_pole,Date)

plot(X_pole);

%selecting the length of the part of the signal will be trasformed
N_ft=N;
%fast Fourier transform calculation
Ftrns_X=fft(X_pole,N_ft);


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

%amplitude spectrum calculation

ampl_spectrum_X=abs(Ftrns_X)/N_ft;

% we plot spectrum multiplied by two to take into account energy in the
% second part
plot(freq, 2*ampl_spectrum_X(2:N_half+1));
title('amplitude spectrum - module of Fourier-transformation ')
xlabel('frequency')


% periods=1./freq;
% plot(periods, 2*apl_spectrum_X(2:N_half+1));
% title('periodogramm')
% xlabel('frequency')
%

%Co L.V. Zotov