Фурье. Загрузить файл с данными

%%
function Lab_Matlab(  )
 
%var
  b=0;
  a=0;
  fi=0;
  tn=0;
 fn=0;
  fd=0;
  td=0;
 sum_ampl=0;
 sum_cos=0;
  sum_sin=0;
sum_cos_fi=0;
 sum_sin_fi=0;
 i=0;
j=0;
n=0;
  ampl=0;
 f=0;
  x=0;
  y=0;
  delta=0;
  z=0;
  sum_x=0;
  Specl=0;
 
b(1)=0.000881847553;
b(2)=-0.00000187613842;
b(3)=0.0127260558;
b(4)=0.00000117333529;
b(5)=-0.0385919636;
b(6)=-0.000000794669754;
b(7)=0.0217106507;
b(8)=-0.000000460638630;
b(9)=0.0994375191;
b(10)=-0.000000126607506;
b(11)=-0.279463243;
b(12)=-0.000000218474142;
b(13)=0.368401231;
b(14)=0.368401231;
b(15)=-0.000000218474142;
b(16)=-0.279463243;
b(17)=-0.000000126607506;
b(18)=0.0994375191;
b(19)=-0.000000460638630;
b(20)=0.0217106507;
b(21)=-0.000000794669754;
b(22)=-0.0385919636;
b(23)=0.00000117333529;
b(24)=0.0127260558;
b(25)=-0.00000187613842;
b(26)=0.000881847553;
fn=0.005;
 
% АЧХ разрабатываемого КИХ-фильтра
i=1;
while i<99 
    sum_cos=0.0;
    sum_sin=0.0;
    j=1;
       while j<26
          sum_cos=sum_cos+b(j)*cos(2*3.14159265358979*fn*j);
          sum_sin=sum_sin+b(j)*sin(2*3.14159265358979*fn*j);
          j=j+1;
        end;
a(i)=sqrt((sum_cos * sum_cos)+(sum_sin * sum_sin));
fn= fn + 0.005;
i=i+1;
end;
str='АЧХ разрабатываемого КИХ-фильтра';
plot(a,'r');
i=1;
 
% ФЧХ разрабатываемого КИХ-фильтра  
fn=0.005;
while i<100 
    sum_cos_fi=0.0;
    sum_sin_fi=0.0;
    j=1;
    while j<26  
        sum_sin_fi=sum_sin_fi+b(j)*sin(2*3.14159265358979*fn*j);
        sum_cos_fi=sum_cos_fi+b(j)*cos(2*3.14159265358979*fn*j);
        j=j+1;
    end;
fi(i)=atan(sum_sin_fi/sum_cos_fi);
fn_list(i)=fn;
fn= fn + 0.005;
i=i+1;
end;
 
figure(2); 
plot(fn_list, fi,'r');
title('ФЧХ разрабатываемого КИХ-фильтра');
 
n=1;
ampl(1)=9;
ampl(2)=11;
ampl(3)=15;
ampl(4)=17;
ampl(5)=21;
ampl(6)=29;
ampl(7)=35;
 
f(1)=21000;
f(2)=25000;
f(3)=32000;
f(4)=33000;
f(5)=41000;
f(6)=45000;
f(7)=50000;
i=1;
fd=2*f(7);
td=1.0/fd;
while n<1024 
    tn(n)=td*n;
    n=n+1;
end;
n=1;
 
% Исходный сигнал 
while n<1024 
    sum_ampl=0.0;
    i=1;
    while i<7
        sum_ampl=sum_ampl+ampl(i)*sin(2*3.14159265358979*f(i)*tn(n));
        i=i+1;
    end;
x(n)=sum_ampl;
n=n+1;
end;
 
figure(3); 
plot(x(1:250),'r');
title('Исходный сигнал');
hold on;
 
% Исходный сигнал с использованием фильтра
z=1;
n=1;
while n<1024 
    sum_x=0.0;
    j=1;
    while j<26
        if j>n 
            z= 1;
        else z= n - j+1;
        end;
        sum_x= sum_x + b(j) * x(z);
        j= j + 1;
    end;
y(n)= sum_x;
n= n + 1;
end;
 
%figure(4);
plot(y(1:240),'b');
title('Исходный сигнал с использованием фильтра');
hold off;
 
% Шумы
n=1;
 
while n<1024 
    delta(n)=x(n)-y(n);
    n=n+1;
end;
figure(5);    
plot(delta(1:250),'r');
title('Шумы');
 
load('X.dat');
Specl=FfT(X);
plot(Specl(1:600),'r');
end