有点忘了下面这段代码是干啥的了,但是和离散傅里叶变换有关,结果和Origion一样:
function [e_i_phi,out1,out2]=DFT_AMP(in1,in2)
%和origion是一样的结果
%e_i_phi就是exp(1j*phi);
%输入in1,时域x轴数据
%输入in2,时域y轴数据
%输出out1,频域x轴
%输出out2,频域y轴(Amplitude)
ax=in1;
ay=in2;
N=length(ax);%向量长度
fs=1/(ax(2)-ax(1));%采样率,相邻两个时间点间隔的倒数
f=(0:N-1)/N;%归一化,步长1/N;
AY=fft(ay,N);%快速傅里叶变换,得到一系列复数AY
phi=unwrap(angle(AY));%unwrap函数可调参数
ff=@(x)exp(1j*x);
e_i_phi=arrayfun(ff,phi);
%0频幅需要除以N,其余频幅除以N/2!!!!
%而真实信号的幅度被平均分到前后各一半所以要乘以2
%除以N以抵消采样数对结果的影响
x=f*fs;%只要时域数据长度一样,采样率一样,则频域的频率坐标一样
y=abs(AY)*2/N;
%%%%%判断N的奇偶%%%%%%%%
if mod(N,2)==0
N_new=(N+2)/2;%偶数
else
N_new=(N+1)/2;%奇数
end
%%%%%%%%%%%%%%%%%%%
y(1)=abs(AY(1))/N;%x和origin对得上,y是origin中的Amplitude
y(N_new)=abs(AY(N_new))/N;
out1=x(1:N_new);out1_1=x(N_new+1:end);
out2=y(1:N_new);out2_1=y(N_new+1:end);
%%%%绘图,可以注释掉%%%%%%%
plot(x(1:N_new),y(1:N_new),'linewidth',1);
xlabel('freq')
ylabel('Amplitude')
%%%%%%%%%%%%%%%%%%%%%
end