MATLAB信号频谱分析.docx
《MATLAB信号频谱分析.docx》由会员分享，可在线阅读，更多相关《MATLAB信号频谱分析.docx（27页珍藏版）》请在冰豆网上搜索。
MATLAB信号频谱分析
题1－周期信号频谱分析
题目：
给定输入信号x（n）为一周期信号，含基波及若干次谐波分量。
要求：
1、给出信号基波频率及谐波个数
2、基波及各个谐波信号的相位和幅值
clear
clc
A=[0.250.50.25];
f=[124];
fi=[00.5*pi0.75*pi];
N=59;
fs=10;
t=（1:
N）./fs;
x_n=A
（1）.*sin（2*pi*f
（1）*t+fi
（1））+A
（2）.*sin（2*pi.*f
（2）*t+fi
（2））+A（3）.*sin（2pi.*f（3）*t+fi（3））;
save（'Signal.mat','x_n'）
程序如下
fs=10;
N=60;
n=1:
N;
t=n/fs;
x=0.25*sin（2*pi*t）+0.5*sin（2*pi*2*t）+0.25*sin（2*pi*4*t）;
figure
（1）
subplot（211）;plot（x）;grid;
y=fft（x,N）;
mag=abs（y）;
k=0:
length（y）-1;
f=fs/N*k;
subplot（212）;
plot（f,mag）;
xlabel（'Frequence（Hz）'）;
ylabel（'Magnitude'）;grid;
title（'对周期信号的一个周期进行频谱分析'）
figure
（2）
x=[xxxxx];
figure
（2）
subplot（211）;plot（x）;grid;
N=length（x）
y=fft（x,N）;
mag=abs（y）;
k=0:
length（y）-1;
f=fs/N*k
subplot（212）;
plot（f,mag）;
xlabel（'Frequence（Hz）'）;
ylabel（'Magnitude'）;grid;
title（'对周期信号的多个周期进行频谱分析'）
题2－数字滤波器设计
题目：
已知输入信号x（t）=0.6sin（200πt）+sin（400πt）+0.3sin（800πt）。
要求:
通过MATLAB编程设计一
1）低通滤波器,输出信号x（t）的最低频率成份.
2）高通滤波器,输出信号x（t）的最高频率成份.
3）带通滤波器,输出信号x（t）的中间频率成份.
4）带阻滤波器,输出信号x（t）的最低+最高频率成份.
（Rp=1dB,Rs=16dB,Fs=3倍最高频率,截止频率、滤波器设计方法及类型自己选择.）
解：
编写程序如下：
Fs=1200;
N=1024;
n=0:
N;
t=n/Fs;
x=0.6*sin（200*pi*t）+sin（400*pi*t）+0.3*sin（800*pi*t）;
y=fft（x,N）;
mag=abs（y）;
k=0:
length（y）-1;
f=Fs/N*k;
figure
（1）
plot（f,mag）;
xlabel（'Frequence（Hz）'）;ylabel（'Magnitude'）;
title（'N=1024'）;grid;
Wp=0.12*pi;Ws=0.2*pi;rp=1;rs=16;Fs=1200;%数字滤波器指标
wp=Wp*Fs;ws=Ws*Fs;%模拟技术指标
[n,wn]=buttord（wp,ws,rp,rs,'s'）;
[z,p,k]=buttap（n）;
[b,a]=zp2tf（z,p,k）;
[b,a]=lp2lp（b,a,wn）;
figure
（2）
h1=freqs（b,a）;
freqs（b,a）;title（'模拟低通的幅频和相频特性图'）;
[bz,az]=impinvar（b,a,Fs）;
figure（3）
freqz（bz,az,256）;title（'数字底通的幅频和相频特性图'）;
y1=filter（bz,az,x）;
y=fft（y1,N）;
mag=abs（y）;
k=0:
length（y）-1;
f=Fs/N*k;
figure（10）
plot（f,mag）;grid;title（'输出信号x（t）的最低频率成份'）
wp=.4*pi;ws=.35*pi;rp=1;rs=16;Fs=1200;%数字滤波器指标
wap=2*Fs*tan（wp/2）;was=2*Fs*tan（ws/2）;
[n,wn]=buttord（wap,was,rp,rs,'s'）;
[z,p,k]=buttap（n）;
[b,a]=zp2tf（z,p,k）;
[bt,at]=lp2hp（b,a,wap）;
w=0:
pi:
4000*pi;
figure（4）
freqs（bt,at,w）;title（'模拟高通的幅频和相频特性图'）;
[bz,az]=bilinear（bt,at,Fs）;
figure（7）
freqz（bz,az,256,1）;title（'数字高通的幅频和相频特性图'）;
y2=filter（bz,az,x）;
y=fft（y2,N）;
mag=abs（y）;
k=0:
length（y）-1;
f=Fs/N*k;
figure（11）
plot（f,mag）;grid;title（'输出信号x（t）的最高频率成份'）
wp=[0.2*pi0.35*pi];ws=[0.12*pi0.4*pi];rp=1;rs=30;Fs=1200;%数子滤波器指标
wap=2*Fs*tan（wp./2）;was=2*Fs*tan（ws./2）;%模拟滤波器指标
[n,wn]=buttord（wap,was,rp,rs,'s'）;
[z,p,k]=buttap（n）;
[bp,ap]=zp2tf（z,p,k）;
bw=wap
（2）-wap
（1）;w0=sqrt（wap
（1）*wap
（2））;
[bs,as]=lp2bp（bp,ap,w0,bw）;
w=0:
pi:
4000*pi;
figure（6）
freqs（bs,as,w）;title（'模拟带通的幅频和相频图'）;
[bz,az]=bilinear（bs,as,Fs）;
figure（7）
freqz（bz,az,256,Fs）;title（'数字带通的幅频和相频图'）
y3=filter（bz,az,x）;
y=fft（y3,N）;
mag=abs（y）;
k=0:
length（y）-1;
f=Fs/N*k;
figure（12）
plot（f,mag）;grid;title（'输出信号x（t）的中间频率成份'）
wp=[0.15*pi0.4*pi];ws=[0.2*pi0.35*pi];rp=1;rs=16;Fs=1200;%数值指标
wap=2*Fs*tan（wp./2）;was=2*Fs*tan（ws./2）;%模拟指标
[n,wn]=buttord（wap,was,rp,rs,'s'）;
[z,p,k]=buttap（n）;
[b,a]=zp2tf（z,p,k）;
bw=wap
（2）-wap
（1）;w0=sqrt（wap
（1）*wap
（2））;
[bt,at]=lp2bs（b,a,w0,bw）;
figure（8）
w=0:
pi:
4000*pi;
freqs（bt,at,w）;title（'模拟带阻的幅频和相频图'）
[bz1,az1]=bilinear（bt,at,Fs）;%模拟到数字
figure（9）
freqz（bz1,az1,256,Fs）;title（'数字带阻的幅频和相频图'）
y4=filter（bz1,az1,x）;
y=fft（y4,N）;
mag=abs（y）;
k=0:
length（y）-1;
f=Fs/N*k;
figure（13）
plot（f,mag）;grid;title（'输出信号x（t）的最低+最高频率成份'）
题3－语音信号采集与处理初步
题目：
1．语音信号的采集
2．语音信号的频谱分析
3．设计数字滤波器和画出频率响应
4．用滤波器对信号进行滤波
5．比较滤波前后语音信号的波形及频谱
6．回放和存储语音信号
解：
编程如下
[x,fs,bits]=wavread（'类似爱情.wav'）;
wavplay（x,fs）;
X=fft（x,4096）;
magX=abs（X）;
angX=angle（X）;
figure
（1）
subplot（221）;plot（x）;title（'原始信号波形'）;
subplot（222）;plot（X）;title（'原始信号频谱'）;
subplot（223）;plot（magX）;title（'原始信号幅值'）;
subplot（224）;plot（angX）;title（'原始信号相位'）;
n=length（x）;
t=0:
1/fs:
（n-1）/fs;
noise=0.5*sin（100*pi*t）+sin（400*pi*t）+0.5*sin（1600*pi*t）;
figure
（2）
plot（t,noise）;grid;
xlabel（'时间（t）'）;ylabel（'幅度（y）'）;title（'噪声的时域图'）;
Noise=fft（noise,n）;
N1=abs（Noise）;
n1=floor（n/2）;
f=（0:
n1）*fs/n;
figure（3）;%画出噪声的频谱图
plot（f,N1（1:
n1+1））;gridon;
xlabel（'频率（f）'）;ylabel（'幅度（N1）'）;title（'噪声的频谱图'）;
x=x（:
1）;
x1=noise+x';
sound（x1,fs）;
wavwrite（x1,16,'类似爱情1.wav'）;
figure（4）%画出加噪声前后的时域比较图
subplot（221）;plot（t,x）;gridon;
xlabel（'时间（t）'）;ylabel（'幅度（x）'）;title（'加噪声前的时域图'）;
subplot（222）;plot（t,x1）;gridon;
xlabel（'时间（t）'）;ylabel（'幅度（Y1）'）;title（'加噪声后的时域图'）;
X=fft（x,n）;%对加噪声前的采样信号1进行傅立叶变换
MagX=abs（X）;
X1=fft（x1,n）;%对加噪声后的采样信号1进行傅立叶变换
MagX1=abs（X1）;
subplot（223）;plot（f,MagX（1:
n1+1））;gridon;
xlabel（'频率（f）'）;ylabel（'幅度（MagX）'）;title（'加噪声前的频谱图'）;
subplot（224）;plot（f,MagX1（1:
n1+1））;gridon;
xlabel（'频率（f）'）;ylabel（'幅度（MagX1）'）;title（'加噪声后的频谱图'）;
%%%%%%%%%%%%%%%构造800Hz的带阻滤波器%%%%%%%%%%%%%
f0_stop1=800;fc=20;Rp=1;Rs=30;
ws_stop1=[（f0_stop1-fc）/（fs/2）（f0_stop1+fc）/（fs/2）];%通带的拐角频率
wp_stop1=[（f0_stop1-5*fc）/（fs/2）（f0_stop1+5*fc）/（fs/2）];%阻带的拐角频率
[N_stop1,wn_stop1]=cheb1ord（wp_stop1,ws_stop1,Rp,Rs）;%求出切比雪夫I型滤波器的阶数N及频率参数wc
[num_stop1,den_stop1]=cheby1（N_stop1,Rp,wn_stop1,'stop'）;%求出切比雪夫I型带阻数字滤波器的传递函数模型系数
[h_stop1,w_stop1]=freqz（num_stop1,den_stop1,fs）%求出离散系统频率响应的数值
figure（5）;
subplot（231）;%画出带阻滤波器的幅频特性图
plot（w_stop1*fs/（2*pi）,20*log10（abs（h_stop1）））;
xlabel（'频率（w_stop1）'）;ylabel（'幅度（h_stop1）'）;title（'带阻滤波器在800Hz处的幅频特性图'）;gridon;
B_stop1=filter（num_stop1,den_stop1,x1）;%对含噪信号x1进行带阻滤波
pause（3）;
sound（B_stop1,fs）;%读出滤去800Hz噪声后的采样信号
wavwrite（B_stop1,fs,16,'类似爱情2.wav'）;
subplot（234）;%画出带阻滤波器滤去800Hz噪声后的时域图
plot（t,B_stop1）;
xlabel（'时间（t））'）;ylabel（'幅度（B_stop1）'）;title（'用带阻滤波器滤去800Hz噪声后的时域图'）;gridon;
%%%%%%%%%%%%%%%构造200Hz的带阻滤波器%%%%%%%%%%%%%%%%%%%%%%
f0_stop2=200;
ws_stop2=[（f0_stop2-fc）/（fs/2）（f0_stop2+fc）/（fs/2）];%通带的拐角频率
wp_stop2=[（f0_stop2-5*fc）/（fs/2）（f0_stop2+5*fc）/（fs/2）];%阻带的拐角频率
[N_stop2,wn_stop2]=cheb1ord（wp_stop2,ws_stop2,Rp,Rs）;%求出切比雪夫I型滤波器的阶数N及频率参数wc
[num_stop2,den_stop2]=cheby1（N_stop2,Rp,wn_stop2,'stop'）;%求出切比雪夫I型带阻数字滤波器的传递函数模型系数
[h_stop2,w_stop2]=freqz（num_stop2,den_stop2,fs）;%求出离散系统频率响应的数值
subplot（232）;
plot（w_stop2*fs/（2*pi）,20*log10（abs（h_stop2）））;
xlabel（'频率（w_stop2）'）;ylabel（'幅度（h_stop2）'）;title（'带阻滤波器在200Hz处的幅频特性图'）;gridon;
B_stop2=filter（num_stop2,den_stop2,B_stop1）;%对含噪信号Y1进行带阻滤波
pause（3）;
sound（B_stop2,fs）;%读出滤去800Hz和200Hz噪声后的采样信号
wavwrite（B_stop2,fs,16,'类似爱情3.wav'）;
subplot（235）;%画出带阻滤波器滤去800Hz和200Hz噪声后的时域图
plot（t,B_stop2）;
xlabel（'时间（t））'）;ylabel（'幅度（B_stop2）'）;title（'带阻滤波器滤去800Hz和200Hz噪声后的时域图'）;gridon;
%%%%%%%%%%%%%%构造50Hz的带阻滤波器%%%%%%%%%%%%%%%%%%%%%%
f0_stop3=50;
ws_stop3=（f0_stop3-fc）/（fs/2）;%阻带的拐角频率
wp_stop3=（f0_stop3+5*fc）/（fs/2）;%通带的拐角频率
Wws=0.0075;
Wwp=0.0375;
[N_stop3,wc_stop3]=cheb1ord（Wwp,Wws,Rp,Rs,'s'）;%求出切比雪夫I型滤波器的阶数N及频率参数wc
[num_stop3,dem_stop3]=cheby1（N_stop3,Rp,wc_stop3,'high'）;%求出切比雪夫I型滤波器高通数字滤波器的传递函数模型系数
[h_stop3,w_stop3]=freqz（num_stop3,dem_stop3,fs）;%求出离散系统频率响应的数值
subplot（233）;%画出高通滤波器的幅频特性图
plot（w_stop3*fs/（2*pi）,（abs（h_stop3）））;
xlabel（'频率（w_stop3）'）;ylabel（'幅度（h_stop3）'）;title（'高通滤波器在50Hz处的幅频特性图'）;gridon;
B_stop3=filter（num_stop3,dem_stop3,B_stop2）;
pause（3）;
sound（B_stop3,fs）;%读出滤去800Hz和200Hz噪声后的采样信号
wavwrite（B_stop3,fs,16,'类似爱情4.wav'）;
subplot（236）;%画出高通滤波器滤去800Hz和200Hz和50Hz噪声后的时域图
plot（t,B_stop3）;
xlabel（'时间（t））'）;ylabel（'幅度（B_stop3）'）;title（'滤去800Hz和200Hz和50Hz噪声后的时域图'）;gridon;
%%%%%%%%%%%%%%%%%%%%%%%%%%用不同的滤波器滤去相应噪声频率后的时域比较%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure（6）;
subplot（221）;plot（t,x1）;gridon;
xlabel（'时间（t）'）;ylabel（'幅度（x1）'）;title（'加噪声后的时域图'）;
subplot（222）;plot（t,B_stop1）;
xlabel（'时间（t））'）;ylabel（'幅度（B_stop1）'）;
title（'用带阻滤波器滤去800Hz噪声后的时域图'）;gridon;
subplot（223）;plot（t,B_stop2）;
xlabel（'时间（t））'）;ylabel（'幅度（B_stop2）'）;
title（'带阻滤波器滤去800Hz和200Hz噪声后的时域图'）;gridon;
subplot（224）;plot（t,B_stop3）;
xlabel（'时间（t））'）;ylabel（'幅度（B_stop3）'）;
title（'高通滤波器滤去800Hz和200Hz和50Hz噪声后的时域图'）;gridon;
%%%%%%%%%%%%%%%%%%%%%%%%%%用不同的滤波器滤去相应噪声频率后的频域比较%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure（7）;
X1=fft（x1,n）;%对加噪声后的采样信号1进行傅立叶变换
MagX1=abs（X1）;
n1=floor（n/2）;
f=（0:
n1）*fs/n;
X1_B_stop1=fft（B_stop1,n）;%对滤去800Hz噪声后的采样信号进行傅立叶变换
MagX1_B_stop1=abs（X1_B_stop1）;
X1_B_stop2=fft（B_stop2,n）;%对滤去800Hz和200Hz噪声后的采样信号进行傅立叶变换
MagX1_B_stop2=abs（X1_B_stop2）;
X1_B_stop3=fft（B_stop3,n）;%对滤去800Hz和200Hz和50Hz噪声后的采样信号进行傅立叶变换
MagX1_B_stop3=abs（X1_B_stop3）;
subplot（221）;plot（f,MagX1（1:
n1+1））;gridon;
xlabel（'频率（f）'）;ylabel（'幅度（X1）'）;title（'加噪声后的频谱图'）;
subplot（2,2,2）;plot（f,MagX1_B_stop1（1:
n1+1））;gridon;
xlabel（'频率（f）'）;ylabel（'幅度（X1_B_stop1）'）;title（'用带阻滤波器滤去800Hz噪声后的的频谱图'）;
subplot（2,2,3）;plot（f,MagX1_B_stop2（1:
n1+1））;gridon;
xlabel（'频率（f）'）;ylabel（'幅度（X1_B_stop2）'）;title（'用带阻滤波器滤去800Hz和200Hz噪声后的的频谱图'）;
subplot（2,2,4）;plot（f,MagX1_B_stop3（1:
n1+1））;gridon;
xlabel（'频率（f）'）;ylabel（'幅度（X1_B_stop3）'）;title（'用高通滤波器滤去800Hz和200Hz和50Hz噪声后的的频谱图'）;
程序二：
[x,fs,bits]=wavread（'toilet.wav'）;
wavplay（x）;
pause（3）;
wavplay（x,11025）;
pause（3）
wavplay（x,44100）;
X=fft（x,4096）;
magX=abs（X）;
angX=angle（X）;
figure
（1）
subplot（221）;plot（x）;title（'原始信号波形'）;
subplot（222）;plot（X）;title（'原始信号频谱'）;
subplot（223）;plot（magX）;title（'原始信号幅值'）;
subplot（224）;plot（angX）;title（'原始信号相位'）;
pause（5）
%N阶高通滤波器
x=wavread（'toilet.wav'）;
sound（x）;
N=5;wc=0.3;
[b,a]=butter（N,wc,'high'）;
X=fft（x）;
figure
（2）
subplot（321）;plot（x）;title（'滤波前信号的波形'）;
subplot（322）;plot（X）;title（'滤波前信号的频谱'）;
y=filter（b,a,x）;
wavwrite（y,'toilet1.wav'）;
wavplay（y）
Y=fft（y）;
subplot（323）;plot（y）;title（'IIR滤波后信号的波形'）;
subplot（324）;plot（Y）;title（'IIR滤波后信号的频谱'）;
z=fftfilt（b,x）;
wavwrite（z,'toilet2.wav'）;
wavplay（z）
Z=fft（z）;
subplot（325）;plot（z）;title（'FIR滤波后信号的波形'）;
subplot（326）;plot（Z）;title（'FIR滤波后信号的频谱'）;
pause（5）
%N阶低通滤波器
x=wavread（'toilet.wav'）;
sound（x）;
X=fft（x）;
figure（3）
subplot（321）;plot（x）;title（'滤波前信号的波形'）;
subplot（322）;plot（X）;title（'滤波前信号的频谱'）
升级会员 