biaozhunfft.m

Matlab实现频谱分析 其幅值和相位已经标准化

clear;clc;close all;
f1=50;
f2=150;
f3=63;
f4=444;
p1=10;
p2=20;
p3=-30;
p4=-40;
fs=1600;
N=1024;
ff=fs/N;
n=0:N-1;
t=n/fs;
s=3*sin(2*pi*f1*t+p1*pi/180)+0.5*sin(2*pi*f2*t+p2*pi/180)+0.3*sin(2*pi*f3*t+p3*pi/180)+0.1*sin(2*pi*f4*t+p4*pi/180);
%load sam s;
w=0.5-0.5*cos(2*pi*n./N);%hanning窗函数
%s=s.*w;%对信号的加窗
plot(t,s);
title('原始信号');
grid on;
figure;%幅度图片
y=fft(s,N);
mag=abs(2*y/N);%换算成实际的幅度
f=0:fs/(N-1):fs;
plot(f(1:N/2),mag(1:N/2));
title('幅度-频率曲线图');
xlabel('频率(Hz)');ylabel('幅值');
grid on;

%寻找谱峰
for i=1:10%10次谐波
k1(i)=floor(f1*i/ff+1);%谱峰附近2点的幅值
k2(i)=k1(i)+1;
if  mag(k1(i))>mag(k2(i))
max(i)=mag(k1(i));
fm(i)=f(k1(i));
end
if mag(k1(i))<mag(k2(i))
max(i)=mag(k2(i));
fm(i)=f(k2(i));
end
end

figure;%相位特性  
for i=1:N/2
 angle(i)=phase(y(i)); %计算相位
 angle(i)=angle(i)*180/pi+90; %换算为角度
end;
plot(f(1:N/2),angle(1:N/2));   %显示相位图
title('相位-频率曲线图');
xlabel('频率(Hz)');ylabel('相位（度）');
grid on;