搜索
%傅里叶变换clc;clearall;closeall;ticFs=128;%采样频率，频谱图的最大频率T=1/Fs;%采样时间，原始信号的时间间隔L=256;%原始信号的长度，即原始离散信号的点数t=(0:L-1)*T;%原始信号的时间取值范围x=7*cos(2*pi*15*t-pi)+3*cos(2*pi*40*t-90*pi/180)+3*cos(2*pi*30*t-90*pi/180);z=7*cos(2*pi*15*t-pi)+3*cos(2*pi*40*t-90*pi/180);z1=6*cos(2*pi*30*t-90*pi/180);z1(1:L/2)=0;z=z+z1;y=x;%+randn(size(t));figure;plot(t,y)title('含噪信号')xlabel('时间（s）')holdonplot(t,z,'r--')N=2^nextpow2(L);%N为使2^N=L的最小幂Y=fft(y,N)/N*2;Z=fft(z,N)/N*2;%快速傅里叶变换之后每个点的幅值是直流信号以外的原始信号幅值的N/2倍（是直流信号的N倍）f=Fs/N*(0:N-1);%频谱图的频率取值范围A=abs(Y);%幅值A1=abs(Z);B=A;%让很小的数置零.B1=A1;A(A10^-10)=0;%A1(A110^-10)=0;P=angle(Y).*A./B;P1=angle(Z).*A1./B1;P=unwrap(P,pi);%初相位值，以除去了振幅为零时的相位值P1=unwrap(P1,pi);figuresubplot(211)plot(f(1:N/2),A(1:N/2))%函数ffs返回值的数据结构具有对称性，因此只取前一半holdonplot(f(1:N/2),A1(1:N/2),'r--')title('幅值频谱')xlabel('频率（HZ）')ylabel('幅值')gridonsubplot(212)stem(f(1:N/2),P(1:N/2))holdonstem(f(1:N/2),P1(1:N/2),'r')title('相位频谱')xlabel('频率（HZ）')ylabel('相位')toc