http://www.ilovematlab.cn/thread-245794-1-1.html
最近在研究matlab画频谱图,在查找资料过程,在matlab中文论坛中看到一篇帖子,感觉比较有用,简单整理了一下。
一、FFT函数
假设原始信号为x,则fft_x=abs(fft(x))得到的是对信号x的进行快速傅立叶变换后的频谱值,只不过这时候横坐标和纵坐标并不是真正的频率值和幅值,横坐标变换成实际频率值的公式为f=(n-1)*Fs/N,其中n为fft_x的横坐标,Fs为采样率,N为x长度(由此可以看出,fft的分辨率实际上就是Fs/N,也就是1/T,它只和采样时间有关),fft_x纵坐标变换成实际信号的幅度值需要做如下变换,y(1)=fft_x(1)/N,y(i)=fft_x(i)*2/N。
二、freqz函数
freqz是求系统的频率响应,调用格式[h,w]=freqz(x),意思是加入一个系统的单位脉冲响应为x,则这个系统的频率响应为[h,w]=freqz(x),h为幅值,w为角频率(0~2pi),实际上就是对x做了一个傅立叶变换。
该频率响应得到的h经过变换后可以得到的值和fft变换纵坐标值相同。而W对应的实际频率要经过变换f=w*fs/2pi。
默认情况下freqz的长度为512。
三、matlab实例及分析
clc;clear;
Fs=256;
N=512;
% Fs=300;
% N=1200; % 若用这套采样频率和采样点数,求出的频率响应就不太对啊
t=0:1/Fs:(N-1)/Fs; %原始信号要从t=0开始向后取值
Signal=(2+3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*75*t+pi*90/180))';
%%%%%以上为产生原始信号程序%%%%
xm=5*cos(2*pi*30*t);
%%%%%画出原始信号波形图%%%%%
subplot(4,1,1);
plot(t,Signal);
title('原始信号');
%axis([0,(N-1)/Fs,-3,5]);
axis([0,0.5,min(Signal),max(Signal)]);
xlabel('t(s)');
ylabel('Signal(V)');
%%%以上为画出原始信号波形图程序%%%
%%%求信号幅度谱与相位谱%%%%
FFT_Signal=abs(fft(Signal)); %单位为 "V"
angle_Signal=angle(fft(Signal)); %单位为 "rad"
Sita=angle_Signal*180/pi; %单位为 "°"
%%%以上为求信号幅度谱与相位谱%%%%
%%%将幅度谱的横坐标改为对应的真实频率值,纵坐标改为对应的真实幅度值%%%%%
i=1:N/2; %matlab中,向量的下标都是从1开始,别和C++里面的搞混了
x=(i-1)*Fs/N; %将横轴改为对应的频率值。
y(i)=FFT_Signal(i)/(N/2);
y(1)=FFT_Signal(1)/N;
subplot(4,1,2);
plot(x,y); %画出幅度谱
title('原始信号FFT,横纵坐标意义已明确');
xlabel('频率(Hz)');
ylabel('幅值(V)');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%画出原始信号相位谱%%%%%%%%%%%%%%%%%%%%
subplot(4,1,4);
i=1:N/2;
f=(i-1)*Fs/N;
for i=1:N/2
if y(i)<=1
half_Sita(i)=0;
else
half_Sita(i)=Sita(i);
end
end
plot(f,half_Sita);
% axis([0,10,min(angle_Signal*180/pi),max(angle_Signal*180/pi)]);
axis([45,80,-180,180]);
title('原始信号相位谱');
xlabel('频率(Hz)');
ylabel('相位(°)');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%看频率响应,即横坐标为0~2pi的情况%%%%%%
[h1,w]=freqz(Signal); %求频率响应,点数自动变为原来的一半,N /2
h2=abs(h1);
i=1:length(h1);
h2(i)=2*h2(i)/length(Signal);
h2(1)=h2(1)/2;
fw=w/pi*Fs/2;
% figure;
subplot(4,1,3);
plot(fw,h2);
axis([0 Fs/2 min(h2) max(h2)]);
xlabel('w(rad)');
ylabel('幅值');
% grid on;
注意:如为除了freqz之外的长度,则freqz在使用时需要指定长度值,
freqz(Signal,1,length(Signal))。
如果直接使用h,w直接画频谱可以直接使用下面语句。
[H,w]=freqz(h);
plot(w*Fs/(2*pi),20*log10(abs(H)));