题1－周期信号频谱分析
题目：给定输入信号x(n)为一周期信号，含基波及若干次谐波分量。
要求：1、给出信号基波频率及谐波个数
2、基波及各个谐波信号的相位和幅值
clear
clc
A=[0.250.50.25];
f=[124];
fi=[00.5pi0.75pi];
N=59;
fs=10;
t=(1:N)./fs;
x_n=A(1).sin(2pif(1)t+fi(1))+A(2).sin(2pi.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.25sin(2pit)+0.5sin(2pi2t)+0.25sin(2pi4t);
figure(1)
subplot(211);plot(x);grid;
y=fft(x,N);
mag=abs(y);
k=0:length(y)-1;
f=fs/Nk;
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/Nk
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.6sin(200pit)+sin(400pit)+0.3sin(800pit);
y=fft(x,N);
mag=abs(y);
k=0:length(y)-1;
f=Fs/Nk;
figure(1)
plot(f,mag);
xlabel(''Frequence(Hz)'');ylabel(''Magnitude'');
title(''N=1024'');grid;
Wp=0.12pi;Ws=0.2pi;rp=1;rs=16;Fs=1200;%数字滤波器指标
wp=WpFs;ws=WsFs;%模拟技术指标
[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/Nk;
figure(10)
plot(f,mag);grid;title(''输出信号x(t)的最低频率成份'')
wp=.4pi;ws=.35pi;rp=1;rs=16;Fs=1200;%数字滤波器指标
wap=2Fstan(wp/2);was=2Fstan(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:4000pi;
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/Nk;
figure(11)
plot(f,mag);grid;title(''输出信号x(t)的最高频率成份'')
wp=[0.2pi0.35pi];ws=[0.12pi0.4pi];rp=1;rs=30;Fs=1200;%数子滤波器指标
wap=2Fstan(wp./2);was=2Fstan(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:4000pi;
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/Nk;
figure(12)
plot(f,mag);grid;title(''输出信号x(t)的中间频率成份'')
wp=[0.15pi0.4pi];ws=[0.2pi0.35pi];rp=1;rs=16;Fs=1200;%数值指标
wap=2Fstan(wp./2);was=2Fstan(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:4000pi;
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/Nk;
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.5sin(100pit)+sin(400pit)+0.5sin(1600pit);
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-5fc)/(fs/2)(f0_stop1+5fc)/(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_stop1fs/(2pi),20log10(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-5fc)/(fs/2)(f0_stop2+5fc)/(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_stop2fs/(2pi),20log10(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+5fc)/(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_stop3fs/(2pi),(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(''滤波前信号的频谱'');
N=5;wc=0.3;
[b,a]=butter(N,wc);
y=filter(b,a,x);
wavwrite(y,''toilet3.wav'');
wavplay(y)
Y=fft(y);
subplot(323);plot(y);title(''IIR滤波后信号的波形'');
subplot(324);plot(Y);title(''IIR滤波后信号的频谱'');
z=fftfilt(b,x);
wavwrite(z,''toilet4.wav'');
wavplay(z)
Z=fft(z);
subplot(325);plot(z);title(''FIR滤波后信号的波形'');
subplot(326);plot(Z);title(''FIR滤波后信号的频谱'');
pause(5)
%2N阶带通滤波器
x=wavread(''toilet.wav'');
sound(x);
N=5;wc=[0.3,0.6];
[b,a]=butter(N,wc);
X=fft(x);
figure(4)
subplot(321);plot(x);title(''滤波前信号的波形'');
subplot(322);plot(X);title(''滤波前信号的频谱'');
y=filter(b,a,x);
wavwrite(y,''toilet5.wav'');
wavplay(y)
Y=fft(y);
subplot(323);plot(y);title(''IIR滤波后信号的波形'');
subplot(324);plot(Y);title(''IIR滤波后信号的频谱'');
z=fftfilt(b,x);
wavwrite(z,''toilet6.wav'');
wavplay(z)
Z=fft(z);
subplot(325);plot(z);title(''FIR滤波后信号的波形'');
subplot(326);plot(Z);title(''FIR滤波后信号的频谱'');
pause(5)
%2N阶带阻滤波器
x=wavread(''toilet.wav'');
sound(x);
N=5;wc=[0.2,0.7];
[b,a]=butter(N,wc,''stop'');
X=fft(x);
figure(5)
subplot(321);plot(x);title(''滤波前信号的波形'');
subplot(322);plot(X);title(''滤波前信号的频谱'');
y=filter(b,a,x);
wavwrite(y,''toilet7.wav'');
wavplay(y)
Y=fft(y);
subplot(323);plot(y);title(''IIR滤波后信号的波形'');
subplot(324);plot(Y);title(''IIR滤波后信号的频谱'');
z=fftfilt(b,x);
wavwrite(z,''toilet8.wav'');
wavplay(z)
Z=fft(z);
subplot(325);plot(z);title(''FIR滤波后信号的波形'');
subplot(326);plot(Z);title(''FIR滤波后信号的频谱'');
|
|
