将序列分解成三角函数并画图

clear
fs=32;%采样率
N=fs5;%采集了3秒的数据
t=0:1/fs:N/fs;%时间序列
%xn=5cos(2pit2)+2cos(2pit3+pi2/2); %函数

load l1.dat;
xn=l1;
N=length(xn);
plot(xn);%显示函数

yn=fft(xn);%傅里叶变换
mag=abs(yn)/(N/2); %傅里叶后的幅值及修正
f=[0:N]*fs/N;%频率序列
plot(f(1:N/2),mag(1:N/2));%显示幅频图
fai=angle(yn);%弧度.可转为度，也可不转
x=0:1/fs:N/fs;
%将幅值分量很小的cos函数不显示
a = 0;
for i = 1:N/2
if mag(i)>= 1 %取幅值大于0.5的频率分量
a = a+1;
mag_NEW(a) = mag(i); %幅值大于0.5的信号分量的幅值
f_NEW(a) = f(i); %幅值大于0.5的信号分量的频率
fai_NEW(a) = fai(i); %幅值大于0.5的信号分量的相位
end
end

for i=1:a
for j=1:length(x)
z(j)=mag_NEW(i)cos(2pi*f_NEW(i)*x(j)+fai_NEW(i));%频率为f(i)的函数曲线
y(j)=f_NEW(i);

end   
figure(3); 
grid on;
plot3(x,y,z);

     text(x(1),y(1),z(1),['(',num2str(x(1)),',',num2str(y(1)),',',num2str(z(1)),')'],'color','r')

hold on;
set(gca,'Xlim',[0,5]);
set(gca,'Ylim',[0,10]);
set(gca,'Zlim',[-5,5]);
xlabel('时间(s)')
ylabel('频率(hz)')
zlabel('幅值(g)')

end