MATLAB作业
方法一:
N=200;dt=0.001;n=1:200;
x=3*sin(2*pi*10*n*dt)+3*sin(2*pi*30*n*dt)+sin(2*pi*40*n*dt)+sin(2*pi*50*n*dt)+sin(2*pi*60*n*dt)+6*sin(2*pi*80*n*dt);%建立时间序列
X=zeros(1,200);                  %给X一个预先的内存空间,提高运行速度
figure(1),plot(n,x);            matlab求傅里叶变换 %画出时间系列图像
%傅里叶变换
for k=1:200                     
    for n=1:200
    X(k)=X(k)+x(n)*exp(-i*2*pi*n*k/N);
    end
end
f=abs(X);                    %对傅里叶变换后的图像取正数部分
figure(2),plot(f);            %画出傅里叶变换后的图像
%滤波
H=ones(200);H(8:14)=0;H(186:193)=0; %建立40、50、60Hz的时间序列所在空间域数值为0的一维矩阵
for k=1:200
    Y(k)=X(k)*H(k);
end
k=1:200;
figure(3);
plot(k,abs(Y)); %画出滤波之后的时间序列
%傅里叶逆变换
y=zeros(1,200)% 给y一个预先的内存空间,提高运行速度
for n=1:200
    for k=1:200
        y(n)=Y(k)*exp(i*2*pi*n*k/N)+y(n);
    end
    y(n)=y(n)/N;
end
figure(4) ;
n=1:200;
plot(n,y,'r-',n,x,'b-');% 画出原时间序列和滤波后的时间序列,可见滤波之后的时间序列振幅明显变小
方法二:
N=200;dt=0.001;n=1:200;
x=3*sin(2*pi*10*n*dt)+3*sin(2*pi*30*n*dt)+sin(2*pi*40*n*dt)+sin(2*pi*50*n*dt)+sin(2*pi*60*n*dt)+6*sin(2*pi*80*n*dt);
X=zeros(1,200);
figure(1),plot(n,x);
X=fft(x,N);
f=abs(X);
figure(2),plot(f);
H=ones(200);H(8:14)=0;H(186:194)=0;
%滤波
for k=1:200
    Y(k)=X(k)*H(k);
end
k=1:200;
figure(3);
plot(k,abs(Y));
%傅里叶逆变换
y=zeros(1,200);
y=ifft(Y,N) ;
figure(4) ;
n=1:200;
plot(n,y,'r-',n,x,'b-');
Figure1:
figure2:
figure3:
Figure4: