matlab频谱分析.docx
《matlab频谱分析.docx》由会员分享，可在线阅读，更多相关《matlab频谱分析.docx（10页珍藏版）》请在冰豆网上搜索。
matlab频谱分析
设计出一套完整的系统，对信号进行频谱分析和滤波处理；
1．产生一个连续信号，包含低频，中频，高频分量，对其进行采样，进行频谱分析，分别设计三种高通，低通，带通滤波器对信号进行滤波处理，观察滤波后信号的频谱。
2.采集一段含有噪音的语音信号（可以录制含有噪音的信号，或者录制语音后再加进噪音信号），对其进行采样和频谱分析，根据分析结果设计出一合适的滤波器滤除噪音信号。
%写上标题
%设计低通滤波器：
[N,Wc]=buttord（）
%估算得到Butterworth低通滤波器的最小阶数N和3dB截止频率Wc
[a,b]=butter（N,Wc）;%设计Butterworth低通滤波器
[h,f]=freqz（）;%求数字低通滤波器的频率响应
figure
（2）;%打开窗口2
subplot（221）;%图形显示分割窗口
plot（f,abs（h））;%绘制Butterworth低通滤波器的幅频响应图
title（巴氏低通滤波器''）;
grid;%绘制带网格的图像
sf=filter（a,b,s）;%叠加函数S经过低通滤波器以后的新函数
subplot（222）;
plot（t,sf）;%绘制叠加函数S经过低通滤波器以后的时域图形
xlabel（'时间（seconds）'）;
ylabel（'时间按幅度'）;
SF=fft（sf,256）;%对叠加函数S经过低通滤波器以后的新函数进行256点的基—2快速傅立叶变换
w=%新信号角频率
subplot（223）;
plot（））;%绘制叠加函数S经过低通滤波器以后的频谱图
title（'低通滤波后的频谱图'）;
%设计高通滤波器
[N,Wc]=buttord（）
%估算得到Butterworth高通滤波器的最小阶数N和3dB截止频率Wc
[a,b]=butter（N,Wc,'high'）;%设计Butterworth高通滤波器
[h,f]=freqz（）;%求数字高通滤波器的频率响应
figure（3）;
subplot（221）;
plot（））;%绘制Butterworth高通滤波器的幅频响应图
title（'巴氏高通滤波器'）;
grid;%绘制带网格的图像
sf=filter（）;%叠加函数S经过高通滤波器以后的新函数
subplot（222）;
plot（t,sf）;;%绘制叠加函数S经过高通滤波器以后的时域图形
xlabel（'Time（seconds）'）;
ylabel（'Timewaveform'）;
w;%新信号角频率
subplot（223）;
plot（））;%绘制叠加函数S经过高通滤波器以后的频谱图
title（'高通滤波后的频谱图'）;
%设计带通滤波器
[N,Wc]=buttord（[）
%估算得到Butterworth带通滤波器的最小阶数N和3dB截止频率Wc
[a,b]=butter（N,Wc）;%设计Butterworth带通滤波器
[h,f]=freqz（）;%求数字带通滤波器的频率响应
figure（4）;
subplot（221）;
plot（f,abs（h））;%绘制Butterworth带通滤波器的幅频响应图
title（'butterbandpassfilter'）;
grid;%绘制带网格的图像
sf=filter（a,b,s）;%叠加函数S经过带通滤波器以后的新函数
subplot（222）;
plot（t,sf）;%绘制叠加函数S经过带通滤波器以后的时域图形
xlabel（'Time（seconds）'）;
ylabel（'Timewaveform'）;
SF=fft（）;%对叠加函数S经过带通滤波器以后的新函数进行256点的基—2快速傅立叶变换
w=（%新信号角频率
subplot（223）;
plot（'））;%绘制叠加函数S经过带通滤波器以后的频谱图
title（'带通滤波后的频谱图'）;
matlab如何做频谱分析
%----------------------------------------------------------------------------------------------------------------
%方案1：
“x=a*cos（2*pi*w*t）”的形式：
%----------------------------------------------------------------------------------------------------------------
%注意：
1.时域的持续时间范围应较大；
% 2.频率w与序列k的对应关系（N为序列总长度）：
w=1/dt*k/N;
% 3.采样频率1/dt应大于w的2倍
% 4.结果曲线的峰值的横坐标对应的就是w和-w值
%----------------------------------------------------------------------------------------------------------------
a=0.75;
w=pi/3;
dt=0.2;
t=[-30*pi:
dt:
30*pi];
N=size（t,2）;
x=a*cos（2*pi*w*t）;
y=fft（x）;
y=fftshift（y）;
figure;
subplot（2,1,1）;
plot（t,x）;
subplot（2,1,2）;
plot（1/dt*（-N/2+1:
N/2）/N,abs（y））;
%----------------------------------------------------------------------------------------------------------------
%方案2：
“x=a*cos（w*t）”的形式：
%----------------------------------------------------------------------------------------------------------------
%注意：
1.时域的持续时间范围应较大；
% 2.频率w与序列k的对应关系（N为序列总长度）：
w=1/dt*2*pi*k/N;
% 3.采样频率1/dt应大于w/（2*pi）的2倍
% 4.结果曲线的峰值的横坐标对应的就是w和-w值
%----------------------------------------------------------------------------------------------------------------
a=0.75;
w=pi/3;
dt=1;
t=[-20*pi:
dt:
20*pi];
N=size（t,2）;
x=a*cos（w*t）;
y=fft（x）;
y=fftshift（y）;
figure;
subplot（2,1,1）;
plot（t,x）;
subplot（2,1,2）;
plot（1/dt*2*pi*（-N/2:
N/2-1）/N,abs（y））;
备注：
由于使用了fftshift，所以得到的频谱序列关于原点对称，如果不需要负半轴的话自行修改一下就ok了
functionf=frequency（x,fs）
dtlen=length（x）;
t=（0:
dtlen-1）/fs;
subplot（211）;
plot（t,x）;
axistight;
y=abs（fft（x））*2/dtlen;
ff=（0:
dtlen/2-1）*fs/dtlen;
subplot（212）;
f=y（1:
floor（dtlen/2））;
plot（ff,f）;
axistight;
return;
Matlab编程实现FFT实践及频谱分析
内容
1．用Matlab产生正弦波,矩形波,以及白噪声信号，并显示各自时域波形图
2．进行FFT变换，显示各自频谱图，其中采样率，频率、数据长度自选
3．做出上述三种信号的均方根图谱,功率图谱,以及对数均方根图谱
4．用IFFT傅立叶反变换恢复信号，并显示恢复的正弦信号时域波形图
源程序
%*************************************************************************%
% FFT实践及频谱分析 %
%*************************************************************************%
%*************************************************************************%
%***************1.正弦波****************%
fs=100;%设定采样频率
N=128;
n=0:
N-1;
t=n/fs;
f0=10;%设定正弦信号频率
%生成正弦信号
x=sin（2*pi*f0*t）;
figure
（1）;
subplot（231）;
plot（t,x）;%作正弦信号的时域波形
xlabel（'t'）;
ylabel（'y'）;
title（'正弦信号y=2*pi*10t时域波形'）;
grid;
%进行FFT变换并做频谱图
y=fft（x,N）;%进行fft变换
mag=abs（y）;%求幅值
f=（0:
length（y）-1）'*fs/length（y）;%进行对应的频率转换
figure
（1）;
subplot（232）;
plot（f,mag）;%做频谱图
axis（[0,100,0,80]）;
xlabel（'频率（Hz）'）;
ylabel（'幅值'）;
title（'正弦信号y=2*pi*10t幅频谱图N=128'）;
grid;
%求均方根谱
sq=abs（y）;
figure
（1）;
subplot（233）;
plot（f,sq）;
xlabel（'频率（Hz）'）;
ylabel（'均方根谱'）;
title（'正弦信号y=2*pi*10t均方根谱'）;
grid;
%求功率谱
power=sq.^2;
figure
（1）;
subplot（234）;
plot（f,power）;
xlabel（'频率（Hz）'）;
ylabel（'功率谱'）;
title（'正弦信号y=2*pi*10t功率谱'）;
grid;
%求对数谱
ln=log（sq）;
figure
（1）;
subplot（235）;
plot（f,ln）;
xlabel（'频率（Hz）'）;
ylabel（'对数谱'）;
title（'正弦信号y=2*pi*10t对数谱'）;
grid;
%用IFFT恢复原始信号
xifft=ifft（y）;
magx=real（xifft）;
ti=[0:
length（xifft）-1]/fs;
figure
（1）;
subplot（236）;
plot（ti,magx）;
xlabel（'t'）;
ylabel（'y'）;
title（'通过IFFT转换的正弦信号波形'）;
grid;
%****************2.矩形波****************%
fs=10;%设定采样频率
t=-5:
0.1:
5;
x=rectpuls（t,2）;
x=x（1:
99）;
figure
（2）;
subplot（231）;
plot（t（1:
99）,x）;%作矩形波的时域波形
xlabel（'t'）;
ylabel（'y'）;
title（'矩形波时域波形'）;
grid;
%进行FFT变换并做频谱图
y=fft（x）;%进行fft变换
mag=abs（y）;%求幅值
f=（0:
length（y）-1）'*fs/length（y）;%进行对应的频率转换
figure
（2）;
subplot（232）;
plot（f,mag）;%做频谱图
xlabel（'频率（Hz）'）;
ylabel（'幅值'）;
title（'矩形波幅频谱图'）;
grid;
%求均方根谱
sq=abs（y）;
figure
（2）;
subplot（233）;
plot（f,sq）;
xlabel（'频率（Hz）'）;
ylabel（'均方根谱'）;
title（'矩形波均方根谱'）;
grid;
%求功率谱
power=sq.^2;
figure
（2）;
subplot（234）;
plot（f,power）;
xlabel（'频率（Hz）'）;
ylabel（'功率谱'）;
title（'矩形波功率谱'）;
grid;
%求对数谱
ln=log（sq）;
figure
（2）;
subplot（235）;
plot（f,ln）;
xlabel（'频率（Hz）'）;
ylabel（'对数谱'）;
title（'矩形波对数谱'）;
grid;
%用IFFT恢复原始信号
xifft=ifft（y）;
magx=real（xifft）;
ti=[0:
length（xifft）-1]/fs;
figure
（2）;
subplot（236）;
plot（ti,magx）;
xlabel（'t'）;
ylabel（'y'）;
title（'通过IFFT转换的矩形波波形'）;
grid;
%****************3.白噪声****************%
fs=10;%设定采样频率
t=-5:
0.1:
5;
x=zeros（1,100）;
x（50）=100000;
figure（3）;
subplot（231）;
plot（t（1:
100）,x）;%作白噪声的时域波形
xlabel（'t'）;
ylabel（'y'）;
title（'白噪声时域波形'）;
grid;
%进行FFT变换并做频谱图
y=fft（x）;%进行fft变换
mag=abs（y）;%求幅值
f=（0:
length（y）-1）'*fs/length（y）;%进行对应的频率转换
figure（3）;
subplot（232）;
plot（f,mag）;%做频谱图
xlabel（'频率（Hz）'）;
ylabel（'幅值'）;
title（'白噪声幅频谱图'）;
grid;
%求均方根谱
sq=abs（y）;
figure（3）;
subplot（233）;
plot（f,sq）;
xlabel（'频率（Hz）'）;
ylabel（'均方根谱'）;
title（'白噪声均方根谱'）;
grid;
%求功率谱
power=sq.^2;
figure（3）;
subplot（234）;
plot（f,power）;
xlabel（'频率（Hz）'）;
ylabel（'功率谱'）;
title（'白噪声功率谱'）;
grid;
%求对数谱
ln=log（sq）;
figure（3）;
subplot（235）;
plot（f,ln）;
xlabel（'频率（Hz）'）;
ylabel（'对数谱'）;
title（'白噪声对数谱'）;
grid;
%用IFFT恢复原始信号
xifft=ifft（y）;
magx=real（xifft）;
ti=[0:
length（xifft）-1]/fs;
figure（3）;
subplot（236）;
plot（ti,magx）;
xlabel（'t'）;
ylabel（'y'）;
title（'通过IFFT转换的白噪声波形'）;
grid;
升级会员 