傅里叶拟合任意一曲线matlab

%  使用说明：用导入数据功能把数据导入，令时间为变量t，位移为变量y
%%data 为自定义数据
close all;
Fs=10;  %采样频率
T=51.2; %采样时间
N=round(T*Fs);  %采样点数
t1=data(1:1:2^9,1);%时间点
y1=data(1:1:2^9,2);%时间点幅值
fake_y=y1-0.05*rand(512,1); %还原信号
err_y=fake_y-y1;
%  plot(t1,y1,t1,fake_y,t1,err_y);

fft_y=fft(y1,N);  %做FFT
Ayy=abs(fft_y);    %求幅值
real_Ayy=Ayy/(N/2);     %真实幅值转换
f=[0:511]*Fs/N;
% 幅度谱
subplot(2,1,1)
plot(f,real_Ayy);
title('y幅值谱');xlabel('频率/Hz');ylabel('幅值/um');grid;
%%
 
a = 0;
for i = 1:N/2
    if real_Ayy(i)>= 0.001 %取幅值大于0.001的频率分量
        a = a+1;
        my_NEW(a) = real_Ayy(i); %幅值大于0.001的信号分量的幅值
        f_NEW(a) = f(i);    %幅值大于0.001的信号分量的频率
        y1_NEW(a) = fft_y(i);    %幅值大于0.001的信号分量的相位
    end
end
phay = angle(y1_NEW);  %求相位并转化为角度(弧度)
%将幅值大于0.001的信号分量的函数表达式写出来
y_NEW = my_NEW(1)*cos(2*pi*f_NEW(1)*t1+phay(1))/2;
for i=2:a
    y_NEW = my_NEW(i)*cos(2*pi*f_NEW(i)*t1+phay(i)) + y_NEW;
end
 
subplot(2,1,2)
plot(t1,y1,'b-', 'LineWidth',4);
hold on
plot(t1,y_NEW, 'r-','LineWidth',2);
legend('原始信号','拟合信号')

%傅里叶频谱图
figure(2)
stem(f(2:250),real_Ayy(2:250),'LineWidth',0.1)
grid on