◆例1设方波的数学模型为：
]5sin 5
13sin 3
1[sin 4)(000 ++
+
=
t t t E
t f T ωωωπ
，基频： T
πω20=
用MATLAB 软件完成该方波的合成设计
◆ MATLAB 源程序
t=-10:0.1:10; %设定一个数组有201个点，方波周期为20 e=5;w=pi/10; %设定方波幅值为5，w 代表w0 m=-5*sign(t); %给定幅值为5的方波函数 y1=(-4*e/pi)*sin(w*t); %计算1次谐波
y3=(-4*e/pi)*(sin(w*t)+sin(3*w*t)/3); %计算3次谐波
y5=(-4*e/pi)*(sin(w*t)+sin(3*w*t)/3+sin(5*w*t)/5); %计算5次谐波 plot(t,y1,'y');hold; grid; %用黄色点线画出1次谐波及网格线,并在同一张图上画其余曲线
plot(t,y3,'g'); %用绿色点线画出3次谐波 plot(t,y5,'b'); %用蓝色点线画出5次谐波 plot(t,m,'-k'); %用黑色实线画方波
title('方波合成');xlabel('t');ylabel('f(t)'); %为图形加上标题 n=50; %合成任意次方波，n 决定方波的合成次数，在此给定50 yn=0; %设置初始值 for i=1:n
yn=yn+(-4*e/pi)*(1/(2*i-1))*sin((2*i-1)*w*t); end; %计算n 次谐波合成
plot(t,yn,'r') %用红色实线画出n 次谐波合成
◆ 从图中我们可以看到Gibbs 现象。在函数的间断点附近，增加傅里叶级数的展开次数，虽然可以使其间断点附近的微小振动的周期变小，但振幅却不能变小。此现象在控制系统表现为：当求控制系统对阶跃函数的响应时，超调量总是存在的。
◆例2（P110）MATLAB中函数FFT应用举例。
%MA TLAB中函数FFT应用举例
t=0:0.001:0.6;
x=sin(2*pi*50*t)+sin(2*pi*120*t);
y=x+2*randn(size(t))
subplot(2,1,1)
plot(y(1:50))
xlable(‘时间轴t’)
ylable(‘信号值f(t)’)
title(‘正弦波＋随即噪声’,’FontSize’,10)
y=fft(y,512);
f=1000*(0:256)/512;
subplot(2,1,2)
plot(f,Y(1:257))
set(gca,’Xtick’,[0,50,100,150,200,250,300,350,400,450,500])
set(gca,’XtickLabel’,’0|50|100|150|200|250|300|350|400|450|500|’) xlabel(‘频率轴\omega’)
ylabel(‘频谱幅值F(\omega)’)
title(‘信号频谱’,’FontSize’,10)
◆例3
例 3.8.3 有一二阶系统，阻尼比ξ＝0.47，固有频率500
ωHz，采样间隔
=
n
=
?t s,采样点数N＝256。试计算理论幅频特性与由系统阶跃响应计算出的幅频特性0004
.0
数据值，并画出两个计算结果的幅频特性曲线。
％example 3.8.3 MA TLAB PROGRAM
N=256;
dt=0.0004
wn=500;
seta=0.47;
dw=2*pi/(N*dt);
a=wn^2;
b=[1,2*seta*wn,a];
t=[0:dt:(N-1)*dt];
c=step(a,b,t);
w=[0:dw:(N-1)*dw];
[mag,phase]=bode(a,b,w);
ycw=fft(c);
Re=real(ycw);
Im=imag(ycw);
for i=1:N
Rw(i,1)=1-Im(i,1)*(i-1)*dw*dt;
Iw(i,1)=Re(i,1)*(i-1)*dw*dt;
end
ffw=Rw+Iw*sqrt(-1);
A w=abs(ffw)
semilogx(w,20*log10(mag),'r-')
axis([100,10000,-30,10])
text(600,12,'幅频特性')
hold on
semilogx(w,20*log10(A w))
axis([100,10000,-30,10])
grid on
◆例4
例6.2.4 用MA TLAB 中的函数XCORR 求出下列两个周期信号的互相关函数，式中的f=10Hz 。
),2sin()(ft t x π= )902sin(5.0)(
+=ft t y π
％例6.2.4 中计算两个周期信号互相关函数的MA TLAB 程序 N=500; Fs=500; n=0:N-1; t=n/Fs;
Lag=200;
x=sin(2*pi*10*t);
y=0.5*sin(2*pi*10*t+90*pi/180); [c,lags]=xcorr(x,y,Lag,'unbiased'); subplot(2,1,1) plot(t,x,'r') hold on
plot(t,y,'b')
xlabel('t');ylabel('x(t)y(t)'); title('原周期信号') grid
hold off
subplot(2,1,2); plot(lags/Fs,c,'k');
xlabel('t');ylabel('Rxy(t)'); title('互相关函数'); grid
◆例5
例 6.2.5 若有信号为
)()2s i n (2)2s i n ()(21t t f t f t x ?ππ++=
式中，,501Hz f =Hz f 1002=，)(t ω为白噪声（用MA TLAB 中的函数产生）。设采样频率2000=Fs ；试用周期法并应用MA TLAB 编程计算，当数据长度分别为2561=N 和10242=N 两种情况下上述信号的功率谱。
%例6.2.5中周期图法计算信号功率谱的MA TLAB 程序 clf
Fs=2000;
%情况1:数据长度N1=256 N1=256; N1fft=256; n1=0:N1-1; t1=n1/Fs; f1=50;
f2=100;
xn1=sin(2*pi*f1*t1)+sin(2*pi*f2*t1)+randn(1,N1); Pxx1=10*log10(abs(fft(xn1,N1fft).^2)/N1); f1=(0:length(Pxx1)-1)*Fs/length(Pxx1); subplot(2,1,1) plot(f1,Pxx1)
ylabel('功率谱(dB)'); title('数据长度N1=256') grid
%情况2:数据长度N2=1024 N2=1024; N2fft=1024; n2=0:N2-1; t2=n2/Fs;
f1=50; f2=100;
xn2=sin(2*pi*f1*t2)+2*sin(2*pi*f2*t2)+randn(1,N2);
Pxx2=10*log10(abs(fft(xn2,N2fft).^2)/N2);
f2=(0:length(Pxx2)-1)*Fs/length(Pxx2);
subplot(2,1,2)
plot(f2,Pxx2)
xlabel('频率(Hz)');
ylabel('功率谱(dB)');
title('数据长度N2=1024')
◆例6（例3.8.1）分别用conv和FFT算法计算序列：x(n)为在区间[0，1]上均匀分布的N 点随机序列，表示为：x(n)＝rand(1,N),h(n)为均值为零、方差为1的N点高斯分布随机序列，表示为：h(n)＝rand(1,L),试求1≤N≤150时的卷积并比较其运算时间。
%例3.8.1 直接卷积和快速卷积的比较
%
conv_time=zeros(1,150);fft_time=zeros(1,150);
%
for N=1:150
tc=0;tf=0; %初始化
L=2*N-1; %加长序列长度
nu=round((log10(L)/log10(2))+0.45);
L=2^nu; %使点数成为2的幂次
for I=1:100
h=randn(1,N); %产生两个随机序列
x=rand(1,N);
t0=clock;
yc=conv(h,x); %直接卷积计算
t1=etime(clock,t0);
tc=tc+t1; %直接卷积运算的时间
t0=clock;
y2=ifft(fft(h,L).*fft(x,L)); %快速卷积计算
t2=etime(clock,t0);
tf=tf+t2; %快速卷积计算的时间
end
%
conv_time(N)=tc/100; %直接卷积计算的平均时间
fft_time(N)=tf/100; %快速卷积计算的平均时间
end
%
n=1:150;subplot(1,1,1); %图形显示上述两种卷积的计算时间
plot(n(25:150),conv_time(25:150),n(25:150),fft_time(25:150))
上述两种卷积的计算时间的比较如图3.8.5所示。
图3.8.5 两种卷积计算时间比较
◆例7（例3.8.2 ）运用FFT 求取矩形脉冲??
?>≤≤=1
,01t 0 ,1)(t t ω 的谱，说明采样频率低
引起的混叠现象。
（1） 先编写有一定通用性的函数文件cftbyfft.m
[cftbyfft.m]
function [AW,f]=cftbyfft(wt,t,flag)
%cftbyfft.m
%本程序采用FFT 计算连续时间Fouie 变换。输出幅频数据对(f,AW)。
%输入量(wt,t)为已经窗口化了的时间函数wt(t)，它们分别是长度为N 的向量。 %对于时限信号，应使该取值时段与窗口长度相比足够小。以提高频率分辨率。 %对于非时限信号，窗口长度的选取应使窗口外的函数值小 %到可忽略，以提高近似精度。
%输入宗量flag 控制输出CFT 的频率范围。 % flag 取非0时(缺省使用)，频率范围在[0,fs); % flag 取0时，频率范围在[-fs/2,fs/2)。 if nargin==2;flag=1;end
N=length(t); %采样点数，应为2的幂次，以求快速。 T=t(length(t))-t(1); %窗口长度 dt=T/N; %时间分辨率 W0=fft(wt); %实施FFT 变换
W=dt*W0; %算得[0,fs]上的N 点CFT 值 df=1/T; %频率分辨率 n=0:1:(N-1)
%把以上计算结果改写到[-fs/2,fs/2]范围
if flag==0
n=-N/2:(N/2-1);
W=fftshift(W); %产生周期序列的频谱
end
f=n*df; %频率分量向量
AW=abs(W); %幅频谱数据向量
if nargout==0
plot(f,AW);grid,xlabel('频率f');ylabel('|w(f)|')
end
（2）运行以下指令，绘制时域波形和幅频谱
M=5; %做2的幂次用。本例把M设得较小，是为了观察混叠。
tend=1; %波形取非零值的时间长度。
T=10; %窗口化长度应足够大，以减小窗口化引起的泄露“旁瓣”效应。N=2^M; %采样点数，取2的幂是为了使FFT运算较快。
dt=T/N; %以上T、N的取值应使N/T＝fs采样频率大于两倍时间波形
%带宽，以克服采样引起的频谱混叠。
%在本例中，据理论分析知W(f=7.5)=Sa(7.5*pi)=1/(7.5
%*pi)<5%.因此，可近似认为本例时间信号带宽为7.5Hz.
n=0:N-1; %采样序列
t=n*dt; %采样点时间序列
w=zeros(size(t,2),1);
Tow=find((tend-t)>0); %产生非零波形时段的相应序列
w(Tow,1)=ones(length(Tow),1); %在窗口时段内定义的完整波形
plot(t,w,'b','LineWidth',2.5),title('Time Waveform');xlabel('t--->')
由上述程序画出的被窗口化的时间波形如图3.8.12所示：
图3.8.12 被窗口化的时间波形
[AW,f]=cftbyfft(w,t,0);
ff=f+eps; %为避免下面指令出现0/0而采取的措施
AWW=abs(sin(pi*ff)./(pi*ff));
plot(f,AW,'b-',ff,AWW,'r:');
title('Aliasing caused by undersampling')
xlabel('f-->');ylabel('|W(f)|'),legend('by FFT','Theoretical')
由上述程序画出的“欠”采样时引起的混叠曲线如图3.8.13所示：
图3.8.13 “欠”采样时引起的混叠
◆例8（例3.8.5）试分别求出含白噪声干扰的正弦信号与白噪声的自相关函数，并对它们的结果进行比较。
%example 3.8.5 MATLAB PROGRAM
N=1000;
n=0:N-1;
Fs=500;
t=n/Fs;
Lag=100; %相关信号的最大延迟量
x1=sin(2*pi*10*t)+0.6*randn(1,length(t)); %含白噪声的正弦信号x1
[c,lags]=xcorr(x1,Lag,'unbiased') %无偏自相关函数的计算
subplot(2,2,1) %画出x1曲线
plot(t,x1)
xlabel('t');ylabel('x1(t)');
title('含白噪声的正弦信号x1');
grid
subplot(2,2,2) %画出x1自相关曲线
plot(lags/Fs,c)
xlabel('t');ylabel('Rxx1(t)');
title('x1的自相关函数');
grid
x2=randn(1,length(t)); %发生白噪声x2
[c,lags]=xcorr(x2,Lag,'unbiased') %白噪声x2的无偏自相关函数
subplot(2,2,3) %画出x2曲线
plot(t,x2)
xlabel('t');ylabel('x2(t)');
title('白噪声x2');
grid
subplot(2,2,4) %画出x2自相关曲线
plot(lags/Fs,c)
xlabel('t');ylabel('Rxx2(t)');
title('x2的自相关函数');
grid
由上述程序画出的信号与自相关函数的曲线如图3.8.19所示：
图3.8.19 信号与自相关函数
由计算结果和和图3.8.19的曲线可以看出：同时具有周期性和白噪声干扰的信号，其自相关函数不但在τ＝0时具有最大值，而且在τ较大时自相关函数仍有明显的周期性，它的频率和原周期信号的周期相同；而无周期的白噪声，自相关函数在τ＝0时也具有最大值，但当τ稍大时迅速衰减至零附近，利用自相关函数的这一特性可用来识别随机信号中是否含有周期成分和它的频率。
◆例9（例4.4.3）试用MA TLAB语言绘制巴特沃思低通模拟滤波器的平方幅频相应曲
线，滤波器的阶数分别为2，5，10，20。
解
％例4.4.3的MA TLAB程序
n=0:0.01:2
for i=1:4
switch i
case 1
N=2;
case 2
N=5;
case 3
N=10;
case 4
N=20;
end
[z,p,k]=buttap(N); %巴特沃思滤波器原型设计函数，n:阶数;z,p,k:
%滤波器零点，极点和增益
[b,a]=zp2tf(z,p,k); % 零极点增益转换为传递函数
[H,w]=freqs(b,a,n); %模拟滤波器频率响应函数
magH2=(abs(H)).^2; %幅度平方函数
hold on
plot(w,magH2);
axis([0201]);
end
xlabel(‘w/wc’) ;
ylabel(‘∣H(jw)∣^2’’);
grid
图4.4.4 不同阶次巴特沃思滤波器的幅度平方函数
◆例10（例4.4.6 ）试用MA TLAB程序，确定一个模拟低通滤波器的阶数N和截止频率
W c。设计指标为：通带边界频率200π，阻带边界频率300π，通带波纹为1dB，在ΩZ 处，幅度衰减大于18 dB。
解
％设计切比雪夫低通滤波器的MA TLAB程序
wp=200*pi;
ws=300*pi;
Rp=1; %通带波纹
Rs=16; %阻带衰减
%计算滤波器阶数
ebs=sqrt(10^(Rp/10)-1); %波纹系数
A=10^(Rs/20); %A为参变量
Wc=wp; %截止频率
Wr=ws/wp;
g=sqrt(A*A-1)/ebs; %g为参变量
N1=log10(g+sqrt(g*g-1))/log10(Wr+sqrt(Wr*Wr-1));
%滤波器阶数计算
N=ceil(N1); %N取整数
运行上述程序后，可得滤波器的截止频率W c和阶数N为
W c=
628.3185
N＝
4
◆例11（例4.4.4）绘制出阶数分别为2，4，6，8的切比雪夫模拟低通滤波器的平方
幅频响应曲线。
解写出其MA TLAB程序如下：
％绘制切比雪夫低通滤波器幅度平方函数的MA TLAB程序
n=0：0.01：2；
for i=1:4
switch i
case 1
N=2；
case 2
N=4；
case 3
N=6；
case 4
N=8；
end
Rp=1; ％滤波器通带波纹系数
[z,p,k]=cheblap(N,Rp); ％设计切比雪夫模拟低通滤波器原型函数，
％z,p,k分别为滤波器的零点、极点和增益[b,a]=zp2tf(z,p,k); ％零点、极点和增益转换为传递函数
[H,w]=freqs(b,a,n); ％模拟滤波器的频率响应
magH2=(abs(H)).^2; ％幅度平方函数
％Output
posplot=[‘10’num2str(i)]; ％定义字符串变量
subplot(posplot)
plot(w,magH2,‘k’);
ylabel(‘|H(jw)^2|’);
title([‘N=’num2str(N)]);
grid
end
其幅度平方函数曲线如图4.4.6所示。
由图4.4.6可以看出：切比雪夫滤波器在通带内具有等波纹起伏特性，而在阻带内则单调下降。阶次越高，特性越接近矩形。
◆例12（例4.4.7）试设计一巴特沃思模拟带通滤波器，设计要求为：通带频率2kHz~3kHz,
两边的过渡带宽为0.5kHz，若通带波纹1dB，阻带衰减大于100dB.
解
％设计巴特沃思带通滤波器MA TLAB程序
％滤波器设计指标
wp=[2000 3000]*2*pi ；
ws=[1500 3500]*2*pi ；
Rp=1;
Rs=100;
%计算阶数和截止频率
[N,Wn]=buttord (wp,ws,Rp,Rs,‵s‵);
Fc=Wn/(2*pi);
％计算滤波器传递函数多项式系数
[b,a]=butter(N,Wn,‵s‵);
%画出滤波器幅频特性
w=linspace(1,4000,1000)*2*pi; %生成线性等间隔向量
H=freqs(b,a,w);
magH=abs(H);
phaH=unwrap(angle(H));
plot(w/(2*pi),20*log10(magH),‵k‵);
xlabel(‵频率（Hz）‵)；
ylabel(‵幅度（dB）‵);
title(‵巴特沃思模拟滤波器‵)
grid on
由上述程序画出的巴特沃思模拟滤波器的幅频特性如图4.4.11所示。
图4.4.11 巴特沃思带通滤波器的幅频特性
◆例13（例4.4.8 ）试设计一个切比雪夫模拟带阻滤波器。要求的指标为：阻带上下边
界频率为2KHz与3KHz，两侧过渡带为0.5KHz，通带波纹1dB,阻带衰减大于60 dB。
解
％设计切比雪夫带阻滤波器MA TLAB程序
％滤波器设计指标
wp=[2000 3000]*2*pi ；
ws=[1500 3500]*2*pi ；
Rp=1;
Rs=60;
%计算阶数和截止频率
[N,Wn]=cheblord (wp,ws,Rp,Rs,‵s ‵); ％计算滤波器传递函数多项式系数 [b,a]=chebyl(N,Rp,Wn,‵stop ‵,‵s ‵); %画出滤波器幅频特性
w=linspace(1,4000,1000)*2*pi; %生成线性等间隔向量 H=freqs(b,a,w); %相位展开 magH=abs(H);
phaH=unwrap(angle(H));
plot(w/(2*pi),20*log10(magH),‵k ‵); xlabel(‵频率（Hz ）‵)； ylabel(‵幅度（dB ）‵); title(‵切比雪夫模拟滤波器‵)
grid on
由上述程序画出的切比雪夫模拟滤波器的幅频特性如图4.4.12所示。
图4.4.12 切比雪夫带通滤波器的幅频特性
◆ 例14（例5.3.4） 一低通滤波器，其通带和阻带的技术指标分别是
dB 25.0R ,2.0p ==πωp dB 50R , 3.0z z ==πω
试用频率抽样法设计一FIR 具有线性相位的滤波器，取 N=20 。
解 由 19, , 1 , 0k , )()()(20
22 ===k
j
d k
N j d e
H e
H k H ππ
当k=2 时，正好是在通带边界频率p ω 处有一个频率抽样点，即
220
22.0?=
=ππωp
下一个抽样点为k=3 ，是阻带上边界频率z ω ，阻带域通带间无过渡带，即 320
23.0?=
=ππωz
则在通带p ωω≤≤0 内有三个抽样点，阻带πωω≤≤p 上共有七个抽样点。从而有 ] 1 , 1 , 0 , , 0 , 1, 1 , 1 [ =k H
由N=20 ，相位常数 , 5.92/)120(2/)1(=-=-=N α 其相位可表示为
??
?
??≤≤-≤≤-=-=19k 10 , )20(95.09k 0 , 95.020
25.9k k k k
πππφ 根据式（5.3.37）得到H(k) ，利用离散傅里叶变换，求得h(n) ，并有频响内插公式最后可得FIR 滤波器的频响)(ωj e H 。其MA TLAB 程序如下：
%例5.3.4中设计FIR 滤波器MA TLAB 程序
N=20；; 1*N)/pi *2( wl ; 1N : 01 ; 2/)1(N alpha =-=-=
；
1N : 1)2/)1floor((N 2k ; )2/)1floor((N : 0kl ;
] 1 , 25.0 , 25.0 , 0 [wdl ; ] 0 , 0 , 1 , 1[Hdr ];
1 , 1, ,15) 1zeros(, 1, 1 , 1 [Hrs -+-=-====
angH = [ -alpha * ( 2 * pi ) / N * kl ,alpha * (2 * pi ) / N * ( N-k2 ) ]; H = Hrs. * exp (j * angH) ;
= real (ifft (H , N )) ;
[ db , mag , pha ,grd , w ] = freqz – m ( h , 1 ) ; subplot (1 , 1 , 1 )
subplot ( 2 , 2 , 1 ) ; plot (wl ( 1 : 11 ) , ‘o ’ , wdl , Hdr ) ;
axis ( [ 0 , 1 , - 0.1 , 1.1 ] ) ; title ( ‘频率样本： N ＝20 ’) xlabel (‘频率（单位：pi ）’) ;ylabel (‘Hr (k)’)
set (gca ,‘XtickMode ’,‘manual ’,‘XTick ’,[0 , 0.2 , 0.3 , 1 ]); set (gca ,‘YtickMode ’,‘manual ’,‘Ytick ’,[0 , 1]); gird subplot(2 , 2 , 2) ; stem(l , h) ; axis([-1 , N ,-0.1 , 0.3])
title (‘单位抽样响应’ ) ; ylabel(‘h(n)’) ; text (N+1 , -0.1 , ‘n ’)
subplot (2 , 2 , 3 ) ; plot (ww/pi , Hr , wl (1:11)/pi , Hrs (1:11)/pi , Hrs(1:11) , ‘o ’); axis ([ 0 , 1 , -0.2 , 1.2 ]) ; title (‘振幅响应’) xlabel (‘频率（单位: pi ）’) ; ylable (‘Hr(w)’)
set (gca , ‘XTickMode ’ , ‘Xtick ’ , [0 , 0.2 , 0.3 , 1])
set (gca ,‘YtickMode ’ , ‘manual ’ , ‘Ytick ’ , [0 , 1]) ; gird subplot (2, 2, 4 ) ; plot (w / pi , db) ; axis ([0 , 1 –200 , 10]) ; gird title (‘幅度响应’) ; xlabel (‘频率（单位：pi ）’) ; ylable (‘分贝数’ ) ;
set (gca , ‘XtickMode ’,‘Manual ’,‘ XTick ’, [0 ; 0.2 ; 0.3 ; 1]) ; set (gca , ‘YtickMode ’,‘Manual ’,‘YTick ’, [-16 ; 0]) ;
set (gca , ‘YtickLableMode ’,‘Manual ’,‘YTickLabels ’, [‘16’ ;’ 0’]) ;
％
上面程序中的fteqz_m 和Hr_Type2 为扩展函数，不在MA TLAB 所带的工具箱内。需要另行编制，函数fteqz_m 为
function [db , mag , pha , grd , w ] = freqz_m (b , a);
%――――――――――――――――――――――――――― %[db , mag , pha , grd , w ] = freqz_m (b , a); % db : 0-π (dB)
% mag : 0-π
% pha : 0-π
% grd : 0-π
% w : 0-π
% b :
% a :
%
[H , w] = freqz (b , a , 1000 ,’ whole’) ;
H = ( H (1 : 1 : 501 ) )’ ;
ｍag = abs ( H ) ;
db = 20 * log 10 ( ( mag + eps ) / max ( mag ) ) ；
pha = angle ( H ) ;
grd = grpdelay (b , a , w ) ;
扩展函数Hr_Type2 的MA TLAB 程序为
%―――――――――――――――――――――――――――
% [Hr , w , b , L ] = Hr_Type2 (h )
% Hr :
% w : 0-
% b :
L : Hr
% h :
%
M = length ( h ) ;
L = M / 2 ;
ｂ= 2 * [h (L :-1 : 1)] ;
n = [1 : 1 : L ] ; n = n – 0.5 ;
w = [0 : 1 : 500]’* pi / 500 ;
H = cos ( w * n )* b’ ;
相应的单位抽样响应h (n)，频响特性)
H（分别用相对幅度和dB 表示）如图5.3.13
(ωj e
所示。
参考文献：
⒈周敏浩编著.信号处理技术基础.北京航空航天大学出版社，2001.
⒉
⒊
⒋
做脑电波信号处理滴嘿嘿。。Matlab addicted Codes %FEATURE EXTRACTER function [features] = EEGfeaturetrainmod(filename,m) a = 4; b = 7; d = 12; e = 30; signals = 0; for index = 1:9; % read in the first ten EEG data because the files are numbered as ha11test01 rather than ha11test1. s = [filename '0' num2str(index) '.dat']; signal = tread_wfdb(s); if signals == 0; signals = signal; else signals = [signals signal]; end end for index = 10:1:m/2; % read in the rest of the EEG training data s = [filename num2str(index) '.dat']; signal = tread_wfdb(s); if signals == 0;
signals = signal; else signals = [signals signal]; end end %%%%% modification just for varying the training testing ratio ------ for index = 25:1:25+m/2; % read in the rest of the EEG training data s = [filename num2str(index) '.dat']; signal = tread_wfdb(s); if signals == 0; signals = signal; else signals = [signals signal]; end end %%%%%end of modification just for varying the training testing ratio----- for l = 1:m % exrating features (power of each kind of EEG wave forms) [Pxx,f]=pwelch(signals(:,l)-mean(signals(:,l)), [], [], [], 200); % relative power fdelta(l) = sum(Pxx(find(f
安装好MATLAB 2012后再安装目录下点击setup.exe 会出现 "查找安装程序类时出错，查找类时出现异常"的错误提示。该错误的解决方法是进入安装目录下的bin 文件夹双击matlab.exe 对安装程序进行激活。这是可以对该matlab.exe 创建桌面快捷方式，以后运行程序是直接双击该快捷方式即可。 信号运算 1、 信号加 MATLAB 实现: x=x1+x2 2、 信号延迟 y(n)=x(n-k) 3、 信号乘 x=x1.*x2 4、 信号变化幅度 y=k*x 5、 信号翻转 y=fliplr(x) 6、 信号采样和 数学描述：y=∑=2 1)(n n n n x MATLAB 实现： y=sum(x(n1:n2)) 7、 信号采样积 数学描述：∏==2 1)(n n n n x y MATLAB 实现： y=prod(x(n1:n2)) 8、 信号能量 数学描述：∑∞ -∞ == n x n x E 2 | )(| MATLAB 实现：Ex=sum(abs(x)^2)
9、 信号功率 数学描述：∑-== 1 2 | )(|1 P N n x n x N MATLAB 实现：Px=sum((abs(x)^2)/N MATLAB 窗函数 矩形窗 w=boxcar(n) 巴特利特窗 w=bartlett(n) 三角窗 w=triang(n) 布莱克曼窗 w=blackman(n) w=blackman(n,sflag) 海明窗 w=haiming(n) W=haiming(n,sflag) sflag 用来控制窗函数首尾的两个元素值，其取值为symmetric 、periodic 汉宁窗 w=hanning(n) 凯塞窗 w=Kaiser(n,beta) ,beta 用于控制旁瓣的高度。n 一定时，beta 越大，其频谱的旁瓣越小，但主瓣宽度相应增加；当beta 一定时，n 发生变化，其旁瓣高度不变。 切比雪夫窗：主瓣宽度最小，具有等波纹型，切比雪夫窗在边沿的采样点有尖峰。 W=chebwin(n,r)
哈尔滨工程大学 实验报告 实验名称：用matlab描述基本信号 班级：电子信息工程4班 学号： 姓名： 实验时间：2016年10月10日 成绩：________________________________ 指导教师：栾晓明 实验室名称：数字信号处理实验室哈尔滨工程大学实验室与资产管理处制
实验一用matlab 描述基本信号 一、 冲激信号 1、 原理： 最简单的信号是（移位的）单位冲激信号： δ[n -n 0] = ? ??≠=00 0 1n n n n （3.1） 在MA TLAB 中产生冲激信号，必须先确定所关注信号部分的长度。如果准备用冲激信 号δ[n ]来激励因果LTI 系统，可能需要观察从n = 0到n = L -1总共L 个点。若选择L = 31，下面的MA TLAB 代码将产生一个“冲激信号”。 1. L = 31; 2. nn = 0 : (L-1); 3. imp = zeros(L, 1); 4. imp(1) = 1; 注意，根据MA TLAB 编址约定，n =0标号必须对应imp(1)。 例：产生移位冲激信号程序（函数文件） function [x,n] = impseq(n0,n1,n2) % 产生 x(n) = delta(n-n0); n1 <=n0 <= n2 % ---------------------------------------------- % [x,n] = impseq(n0,n1,n2) % if ((n0 < n1) | (n0 > n2) | (n1 > n2)) error('参数必须满足 n1 <= n0 <= n2') end n = [n1:n2]; %x = [zeros(1,(n0-n1)), 1, zeros(1,(n2-n0))]; x = [(n-n0) == 0]; 以上函数文件可以产生指定区间内的冲激移位脉冲。 例1—1：调用这个函数文件生成并绘制： x(n) = 2δ[n+2]-δ[n -4] -5≤ n ≤ 5 程序 % x(n) = 2*delta(n+2) - delta(n-4), -5<=n<=5 n = [-5:5]; x = 2*impseq(-2,-5,5)-impseq(4,-5,5); stem(n,x); title('例 2.1a 的序列图') ylabel('x(n)'); axis([-5,5,-2,3]);text(5.5,-2,'n')
◆例1设方波的数学模型为： ]5sin 513sin 31[sin 4)(000 t t t E t f T ，基频： T 20 用MATLAB 软件完成该方波的合成设计 ◆ MATLAB 源程序 t=-10:0.1:10; %设定一个数组有201个点，方波周期为20 e=5;w=pi/10; %设定方波幅值为5，w 代表w0 m=-5*sign(t); %给定幅值为5的方波函数 y1=(-4*e/pi)*sin(w*t); %计算1次谐波 y3=(-4*e/pi)*(sin(w*t)+sin(3*w*t)/3); %计算3次谐波 y5=(-4*e/pi)*(sin(w*t)+sin(3*w*t)/3+sin(5*w*t)/5); %计算5次谐波 plot(t,y1,'y');hold; grid; %用黄色点线画出1次谐波及网格线,并在同一张图上画其余曲线 plot(t,y3,'g'); %用绿色点线画出3次谐波 plot(t,y5,'b'); %用蓝色点线画出5次谐波 plot(t,m,'-k'); %用黑色实线画方波 title('方波合成');xlabel('t');ylabel('f(t)'); %为图形加上标题 n=50; %合成任意次方波，n 决定方波的合成次数，在此给定50 yn=0; %设置初始值 for i=1:n yn=yn+(-4*e/pi)*(1/(2*i-1))*sin((2*i-1)*w*t); end; %计算n 次谐波合成 plot(t,yn,'r') %用红色实线画出n 次谐波合成 ◆ 从图中我们可以看到Gibbs 现象。在函数的间断点附近，增加傅里叶级数的展开次数，虽然可以使其间断点附近的微小振动的周期变小，但振幅却不能变小。此现象在控制系统表现为：当求控制系统对阶跃函数的响应时，超调量总是存在的。
语音信号处理系统设计 摘要：语音信号处理是研究用数字信号处理技术对语音信号进行处理的一门学科。语音信号处理的目的是得到某些参数以便高效传输或存储,或者是用于某种应用，如人工合成出语音、辨识出讲话者、识别出讲话内容、进行语音增强等。本文简要介绍了语音信号采集与分析以及语音信号的特征、采集与分析方法，并在采集语音信号后，在MATLAB 软件平台上进行频谱分析,并对所采集的语音信号加入干扰噪声，对加入噪声的信号进行频谱分析，设计合适的滤波器滤除噪声，恢复原信号。利用MATLAB来读入（采集）语音信号，将它赋值给某一向量，再将该向量看作一个普通的信号，对其进行FFT变换实现频谱分析，再依据实际情况对它进行滤波，然后我们还可以通过sound命令来对语音信号进行回放，以便在听觉上来感受声音的变化。 关键词：Matlab，语音信号，傅里叶变换，滤波器 1课程设计的目的和意义 本设计课题主要研究语音信号初步分析的软件实现方法、滤波器的设计及应用。通过完成本课题的设计，拟主要达到以下几个目的： 1.1．了解Matlab软件的特点和使用方法。 1.2．掌握利用Matlab分析信号和系统的时域、频域特性的方法； 1.3．掌握数字滤波器的设计方法及应用。 1.4．了解语音信号的特性及分析方法。 1.5．通过本课题的设计，培养学生运用所学知识分析和解决实际问题的能力。 2 设计任务及技术指标 设计一个简单的语音信号分析系统，实现对语音信号时域波形显示、进行频谱分析，
利用滤波器滤除噪声、对语音信号的参数进行提取分析等功能。采用Matlab设计语言信号分析相关程序，并且利用GUI设计图形用户界面。具体任务是： 2.1．采集语音信号。 2.2．对原始语音信号加入干扰噪声，对原始语音信号及带噪语音信号进行时频域分析。 2.3．针对语音信号频谱及噪声频率，设计合适的数字滤波器滤除噪声。 2.4．对噪声滤除前后的语音进行时频域分析。 2.5.对语音信号进行重采样，回放并与原始信号进行比较。 2.6．对语音信号部分时域参数进行提取。 2.7．设计图形用户界面（包含以上功能）。 3 设计方案论证 3.1语音信号的采集 使用电脑的声卡设备采集一段语音信号，并将其保存在电脑中。 3.2语音信号的处理 语音信号的处理主要包括信号的提取播放、信号的重采样、信号加入噪声、信号的傅里叶变换和滤波等，以及GUI图形用户界面设计。 Ⅰ.语音信号的时域分析 语音信号是一种非平稳的时变信号，它携带着各种信息。在语音编码、语音合成、语音识别和语音增强等语音处理中无一例外需要提取语音中包含的各种信息。语音信号分析的目的就在与方便有效的提取并表示语音信号所携带的信息。语音信号分析可以分为时域和变换域等处理方法，其中时域分析是最简单的方法。 Ⅱ.语音信号的频域分析 信号的傅立叶表示在信号的分析与处理中起着重要的作用。因为对于线性系统来说，可以很方便地确定其对正弦或复指数和的响应，所以傅立叶分析方法能完善地解决许多信号分析和处理问题。另外，傅立叶表示使信号的某些特性变得更明显，因此，它能更
系（院）物理与电子工程学院专业电子信息工程题目语音信号的处理与分析 学生姓名 指导教师 班级 学号 完成日期：2013 年5 月 目录 1 绪论 (3) 1.1课题背景及意义 (3) 1.2国内外研究现状 (3) 1.3本课题的研究内容和方法 (4) 1.3.1 研究内容 (4) 1.3.2 开发环境 (4) 2 语音信号处理的总体方案 (4) 2.1 系统基本概述 (4) 2.2 系统基本要求与目的 (4) 2.3 系统框架及实现 (5) 2.3.1 语音信号的采样 (5) 2.3.2 语音信号的频谱分析 (5) 2.3.3 音乐信号的抽取 (5) 2.3.4 音乐信号的AM调制 (5) 2.3.5 AM调制音乐信号的同步解调 (5) 2.4系统设计流程图 (6) 3 语音信号处理基本知识 (6) 3.1语音的录入与打开 (6)
3.2采样位数和采样频率 (6) 3.3时域信号的FFT分析 (7) 3.4切比雪夫滤波器 (7) 3.5数字滤波器设计原理 (8) 4 语音信号实例处理设计 (8) 4.1语音信号的采集 (8) 4.3.1高频调制与低频调制 (10) 4.3.2切比雪夫滤波 (11) 4.3.3 FIR滤波 (11) 5 总结 (12) 参考文献 (13) 语音信号的处理与分析 【摘要】语音信号处理是研究用数字信号处理技术和语音学知识对语音信号进行处理的新兴的学科，是目前发展最为迅速的信息科学研究领域的核心技术之一。通过语音传递信息是人类最重要、最有效、最常用和最方便的交换信息形式。 Matlab语言是一种数据分析和处理功能十分强大的计算机应用软件，它可以将声音文件变换为离散的数据文件，然后利用其强大的矩阵运算能力处理数据，如数字滤波、傅里叶变换、时域和频域分析、声音回放以及各种图的呈现等，它的信号处理与分析工具箱为语音信号分析提供了十分丰富的功能函数，利用这些功能函数可以快捷而又方便地完成语音信号的处理和分析以及信号的可视化，使人机交互更加便捷。信号处理是Matlab重要应用的领域之一。 本设计针对现在大部分语音处理软件内容繁多、操作不便等问题，采用MATLAB7.0综合运用GUI界面设计、各种函数调用等来实现语音信号的变频、变幅、傅里叶变换及滤波，程序界面简练，操作简便，具有一定的实际应用意义。 最后，本文对语音信号处理的进一步发展方向提出了自己的看法。 【关键词】Matlab 语音信号傅里叶变换低通滤波器
MATLAB 在数字信号处理中的应用：连续信号的采样与重建 一、 设计目的和意义 随着通信技术的迅速发展以及计算机的广泛应用，利用数字系统处理模拟信号的情况变得更加普遍。数字电子计算机所处理和传送的都是不连续的数字信号，而实际中遇到的大都是连续变化的模拟量，现代应用中经常要求对模拟信号采样，将其转换为数字信号，然后对其进行计算处理，最好在重建为模拟信号。 采样在连续时间信号与离散时间信号之间其桥梁作用，是模拟信号数字化的第一个步骤，研究的重点是确定合适的采样频率，使得既要能够从采样信号（采样序列）中五失真地恢复原模拟信号，同时由要尽量降低采样频率，减少编码数据速率，有利于数据的存储、处理和传输。 本次设计中，通过使用用MATLAB 对信号f （t ）=A1sin(2πft)+A2sin(4πft)+A3sin(5πft)在300Hz 的频率点上进行采样，并进行仿真，进一步了解MA TLAB 在数字信号处理上的应用，更加深入的了解MA TLAB 的功能。 二、 设计原理 1、 时域抽样定理 令连续信号 xa(t)的傅立叶变换为Xa （j Ω），抽样脉冲序列p(t)傅立叶变换为P （j Ω），抽样后的信号x^(t)的傅立叶变换为X^(j Ω)若采用均匀抽样，抽样周期Ts ，抽样频率为Ωs= 2πfs ，有前面分析可知：抽样过程可以通过抽样脉冲序列p （t ）与连续信号xa （t ）相乘来完成，即满足：x^(t)p(t),又周期信号f （t ）傅立叶变换为： F[f(t)]=2[(]n s n F j n π δ∞ =-∞Ω-Ω∑ 故可以推得p(t)的傅立叶变换为： P （j Ω）=2[(]n s n P j n π δ∞ =-∞Ω-Ω∑ 其中： 根据卷积定理可知： X （j Ω）=12π Xa （j Ω）*P(j Ω) 得到抽样信号x （t ）的傅立叶变换为： X （j Ω）= [()]n n s n P X j n ∞=-∞Ω-Ω∑ 其表明：信号在时域被抽样后，他的频率X （j Ω）是连续信号频率X （j Ω）的形状以抽样频率Ωs 为间隔周期重复而得到，在重复过程中幅度被p （t ）的傅立叶级数Pn 加权。因为只是n 的函数，所以X （j Ω）在重复过程中不会使其形状发生变化。 假定信号x （t ）的频谱限制在-Ωm~+Ωm 的范围内，若以间隔Ts 对xa （t ）进行抽样信号X^（j Ω）是以Ωs 为周期重复。显然，若早抽样过程中Ωs<Ωm ，则 X^ (j Ω)将会发生频谱混叠的现象，只有在抽样的过程中满足Ωs>2Ωm 条件，X^(j Ω)才不会产生混频的混叠，在接收端完全可以有x^（t ）恢复原连续信号xa （t ），这就是低通信号的抽样定理的核心内容。
短时能量分析matlab源程序： x=wavread('4.wav'); %计算N=50，帧移=50时的语音能量 s=fra(50,50,x);%对输入的语音信号进行分帧，其中帧长50，帧移50 s2=s.^2;%一帧内各种点的能量 energy=sum(s2,2);%求一帧能量 subplot(2,2,1); plot(energy) xlabel('帧数'); ylabel('短时能量E'); legend('N=50'); axis([0,500,0,30]) %计算N=100，帧移=100时的语音能量 s=fra(100,100,x); s2=s.^2; energy=sum(s2,2); subplot(2,2,2); plot(energy) xlabel('帧数'); ylabel('短时能量E'); legend('N=100'); axis([0,300,0,30]) %计算N=400，帧移=400时的语音能量 s=fra(400,400,x); s2=s.^2; energy=sum(s2,2); subplot(2,2,3); plot(energy) xlabel('帧数'); ylabel('短时能量E'); legend('N=400'); axis([0,60,0,100]) %计算N=800，帧移=800时的语音能量 s=fra(800,800,x); s2=s.^2; energy=sum(s2,2); subplot(2,2,4); plot(energy) xlabel('帧数'); ylabel('短时能量E'); legend('N=800'); axis([0,30,0,200]) 分帧子函数： function f=fra(len,inc,x) %对读入语音分帧，len为帧长，inc为帧重叠样点数，x为输入语音数据 fh=fix(((size(x,1)-len)/inc)+1);%计算帧数 f=zeros(fh,len);%设一个零矩阵，行为帧数，列为帧长 i=1;n=1; while i<=fh %帧间循环 j=1; while j<=len %帧内循环 f(i,j)=x(n); j=j+1;n=n+1; end n=n-len+inc;%下一帧开始位置 i=i+1; end
工程设计论文 题目：基于MATLAB的语音信号采集与处理 姓名： 班级： 学号： 指导老师：
一.选题背景 1、实践意义： 语音信号是一种非平稳的时变信号，它携带着各种信息。在语音编码、语音合成、语音识别和语音增强等语音处理中无一例外需要提取语音中包含的各种信息。语音信号分析的目的就在于方便有效地提取并表示语音信号所携带的信息。所以理解并掌握语音信号的时域和频域特性是非常重要的。 通过语音相互传递信息是人类最重要的基本功能之一.语言是人类特有的功能.声音是人类常用工具,是相互传递信息的最重要的手段.虽然,人可以通过多种手段获得外界信息,但最重要,最精细的信息源只有语言,图像和文字三种.与用声音传递信息相比,显然用视觉和文字相互传递信息,其效果要差得多.这是因为语音中除包含实际发音容的话言信息外,还包括发音者是谁及喜怒哀乐等各种信息.所以,语音是人类最重要,最有效,最常用和最方便的交换信息的形式.另一方面,语言和语音与人的智力活动密切相关,与文化和社会的进步紧密相连,它具有最大的信息容量和最高的智能水平。 语音信号处理是研究用数字信号处理技术对语音信号进行处理的一门学科,处理的目的是用于得到某些参数以便高效传输或存储;或者是用于某种应用,如人工合成出语音,辨识出讲话者,识别出讲话容,进行语音增强等. 语音信号处理是一门新兴的学科,同时又是综合性的多学科领域,
是一门涉及面很广的交叉学科.虽然从事达一领域研究的人员主要来自信息处理及计算机等学科.但是它与语音学,语言学,声学,认知科学,生理学,心理学及数理统计等许多学科也有非常密切的联系. 语音信号处理是许多信息领域应用的核心技术之一,是目前发展最为迅速的信息科学研究领域中的一个.语音处理是目前极为活跃和热门的研究领域,其研究涉及一系列前沿科研课题,巳处于迅速发展之中;其研究成果具有重要的学术及应用价值. 数字信号处理是利用计算机或专用处理设备，以数值计算的方法对信号进行采集、抽样、变换、综合、估值与识别等加工处理，借以达到提取信息和便于应用的目的。它在语音、雷达、图像、系统控制、通信、航空航天、生物医学等众多领域都获得了极其广泛的应用。具有灵活、精确、抗干扰强、度快等优点。 数字滤波器, 是数字信号处理中及其重要的一部分。随着信息时代和数字技术的发展，受到人们越来越多的重视。数字滤波器可以通过数值运算实现滤波，所以数字滤波器处理精度高、稳定、体积小、重量轻、灵活不存在阻抗匹配问题，可以实现模拟滤波器无法实现的特殊功能。数字滤波器种类很多，根据其实现的网络结构或者其冲激响应函数的时域特性，可分为两种，即有限冲激响应( FIR，Finite Impulse Response)滤波器和无限冲激响应( IIR，Infinite Impulse Response)滤波器。 FIR滤波器结构上主要是非递归结构，没有输出到输入的反馈，系统函数H (z)在处收敛，极点全部在z = 0处（因果系统），因而只能
Matlab滤波信号处理函数 2009-12-04 19:32:22| 分类：matlab方法| 标签：|字号大中小订阅 1 conv 功能：求卷积。 格式：c = conv(a,b) 说明：c = conv(a,b)返回向量a、b的卷积c。 举例：a = [1 2 3] b = [4 5 6] c = conv(a,b) c= 4 13 28 27 18 2 impz 功能：数字滤波器的冲激响应。 格式：[h,t] = impz(b,a) [h,t] = impz(b,a,n) [h,t] = impz(b,a,n,Fs) impz(b,a) impz(...) 说明：[h,t] = impz(b,a)返回系统(b,a)的冲激响应h和相应的时间轴向量t，b、a分别为系统传递函数的分子和分母系数向量。
[h,t] = impz(b,a,n)返回指定的n点冲激响应 [h,t] = impz(b,a,n,Fs)指定了冲激响应采样点的频率间隔1/Fs。Fs 为相对频率， 缺省值为1。 impz(b,a)和impz(...)绘制冲激响应的图形。 举例：计算线性系统(b,a)的冲激响应，结果见图1.4.1。 b =[0.2 0.1 0.3 0.1 0.2]; a =[1 ?.1 1.55 ?.7 0.3]; impz(b,a,50) 3 zplane 功能：离散系统的零极点图。 格式：zplane(z,p) zplane(b,a) 说明：zplane(z,p)和zplane(b,a)绘制系统的零极点图，用“o”表示零点，“x”表示 极点。z、p分别为零点和极点向量，b、a分别为系统传递函数的分子和分母 系数向量。 举例：计算线性系统(b,a)的零点和极点，结果见图1.4.2。 b =[0.2 0.1 0.3 0.1 0.2]; a =[1.0 -1.1 1.5 -0.7 0.3]; zplane(b,a)
MATLAB中FFT的使用方法 说明：以下资源来源于《数字信号处理的MATLAB实现》万永革主编 一.调用方法 X=FFT(x)； X=FFT(x，N)； x=IFFT(X); x=IFFT(X,N) 用MATLAB进行谱分析时注意： （1）函数FFT返回值的数据结构具有对称性。 例： N=8; n=0:N-1; xn=[43267890]; Xk=fft(xn) → Xk= 39.0000-10.7782+6.2929i0-5.0000i 4.7782-7.7071i 5.0000 4.7782+7.7071i0+5.0000i-10.7782-6.2929i Xk与xn的维数相同，共有8个元素。Xk的第一个数对应于直流分量，即频率值为0。 （2）做FFT分析时，幅值大小与FFT选择的点数有关，但不影响分析结果。在IFFT时已经做了处理。要得到真实的振幅值的大小，只要将得到的变换后结果乘以2除以N即可。 二.FFT应用举例 例1：x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)。采样频率fs=100Hz，分别绘制N=128、1024点幅频图。
clf; fs=100;N=128;%采样频率和数据点数 n=0:N-1;t=n/fs;%时间序列 x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);%信号 y=fft(x,N);%对信号进行快速Fourier变换 mag=abs(y);%求得Fourier变换后的振幅 f=n*fs/N;%频率序列 subplot(2,2,1),plot(f,mag);%绘出随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=128');grid on; subplot(2,2,2),plot(f(1:N/2),mag(1:N/2));%绘出Nyquist频率之前随频率变化的振幅xlabel('频率/Hz'); ylabel('振幅');title('N=128');grid on; %对信号采样数据为1024点的处理 fs=100;N=1024;n=0:N-1;t=n/fs; x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);%信号 y=fft(x,N);%对信号进行快速Fourier变换 mag=abs(y);%求取Fourier变换的振幅 f=n*fs/N; subplot(2,2,3),plot(f,mag);%绘出随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=1024');grid on; subplot(2,2,4) plot(f(1:N/2),mag(1:N/2));%绘出Nyquist频率之前随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=1024');grid on; 运行结果：
实验一 基于Matlab 的数字信号处理基本操作 一、 实验目的：学会运用MA TLAB 表示的常用离散时间信号；学会运用MA TLAB 实现离 散时间信号的基本运算。 二、 实验仪器：电脑一台，MATLAB6.5或更高级版本软件一套。 三、 实验内容： （一） 离散时间信号在MATLAB 中的表示 离散时间信号是指在离散时刻才有定义的信号，简称离散信号，或者序列。离散序列通常用)(n x 来表示，自变量必须是整数。 离散时间信号的波形绘制在MATLAB 中一般用stem 函数。stem 函数的基本用法和plot 函数一样，它绘制的波形图的每个样本点上有一个小圆圈，默认是空心的。如果要实心，需使用参数“fill ”、“filled ”，或者参数“.”。由于MATLAB 中矩阵元素的个数有限，所以MA TLAB 只能表示一定时间范围内有限长度的序列；而对于无限序列，也只能在一定时间范围内表示出来。类似于连续时间信号，离散时间信号也有一些典型的离散时间信号。 1. 单位取样序列 单位取样序列)(n δ，也称为单位冲激序列，定义为 ) 0() 0(0 1)(≠=?? ?=n n n δ 要注意，单位冲激序列不是单位冲激函数的简单离散抽样，它在n =0处是取确定的值1。在MATLAB 中，冲激序列可以通过编写以下的impDT .m 文件来实现，即 function y=impDT(n) y=(n==0); %当参数为0时冲激为1,否则为0 调用该函数时n 必须为整数或整数向量。 【实例1-1】 利用MATLAB 的impDT 函数绘出单位冲激序列的波形图。 解：MATLAB 源程序为 >>n=-3:3; >>x=impDT(n); >>stem(n,x,'fill'),xlabel('n'),grid on >>title('单位冲激序列') >>axis([-3 3 -0.1 1.1]) 程序运行结果如图1-1所示。 图1-1 单位冲激序列
郑州轻工业学院 课程设计任务书 题目利用MATLAB实现信号的AM调制与解调专业、班级电子信息工程级班学号姓名 主要内容、基本要求、主要参考资料等： 主要内容： 利用MATLAB对信号 () () ?? ? ? ?≤ = 其他 ,0 t , 100 2t t Sa t m 进行AM调制，载波信号 频率为1000Hz，调制深度为0.5。t0=0.2;首先在MATLAB中显示调制信号的波形和频谱，已调信号的波形和频谱，比较信号调制前后的变化。然后对已调信号解调，并比较解调后的信号与原信号的区别。 基本要求： 1、掌握利用MATLAB实现信号AM调制与解调的方法。 2、学习MATLAB中信号表示的基本方法及绘图函数的调用，实现对常用连续时间信号的可视化表示。 3、加深理解调制信号的变化；验证信号调制的基本概念、基本理论，掌握信号与系统的分析方法。 主要参考资料： 1、王秉钧等. 通信原理[M].北京：清华大学出版社，2006.11 2、陈怀琛.数字信号处理教程----MATLAB释义与实现[M].北京：电子工业出版社,2004. 完成期限：2014.6.9—2014.6.13 指导教师签名： 课程负责人签名： 2014年6月5日
目录 摘要 (1) 1.matlab简介 (2) 1.1matlab基本功能 (2) 1.2matlab应用 (2) 2.系统总体设计方案 (4) 2.1调制信号 (4) 2.1.1 matlab实现调制信号的波形 (4) 2.1.2 matlab实现调制信号的频谱 (4) 2.1.3 matlab实现载波的仿真 (5) 2.2信号的幅度调制 (6) 2.2.1信号的调制 (6) 2.2.2幅度调制原理 (6) 2.2.3 matlab实现双边带幅度调制 (8) 2.2.4 matlab实现已调信号的频谱图 (8) 2.2.5 幅度调制前后的比较 (9) 2.3已调信号的解调 (9) 2.3.1 AM信号的解调原理及方式 (9) 2.3.2 matlab实现已调信号的解调 (11) 2.3.3信号解调前后的比较 (12) 结论与展望 (13) 参考文献 (14) 附录 (15)
实验一 数字信号处理的Matlab 仿真 一、实验目的 1、掌握连续信号及其MA TLAB 实现方法； 2、掌握离散信号及其MA TLAB 实现方法 3、掌握离散信号的基本运算方法，以及MA TLAB 实现 4、了解离散傅里叶变换的MA TLAB 实现 5、了解IIR 数字滤波器设计 6、了解FIR 数字滤波器设计1 二、实验设备 计算机，Matlab 软件 三、实验内容 （一）、 连续信号及其MATLAB 实现 1、 单位冲击信号 ()0,0()1,0 t t t dt εεδδε-?=≠??=?>??? 例1.1：t=1/A=50时，单位脉冲序列的MA TLAB 实现程序如下： clear all; t1=-0.5:0.001:0; A=50; A1=1/A; n1=length(t1); u1=zeros(1,n1); t2=0:0.001:A1; t0=0; u2=A*stepfun(t2,t0); t3=A1:0.001:1; n3=length(t3); u3=zeros(1,n3); t=[t1 t2 t3]; u=[u1 u2 u3]; plot(t,u) axis([-0.5 1 0 A+2]) 2、 任意函数 ()()()f t f t d τδττ+∞ -∞=-? 例1.2：用MA TLAB 画出如下表达式的脉冲序列 ()0.4(2)0.8(1) 1.2() 1.5(1) 1.0(2)0.7(3)f n n n n n n n δδδδδδ=-+-+++++++
clear all; t=-2:1:3; N=length(t); x=zeros(1,N); x(1)=0.4; x(2)=0.8 x(3)=1.2; x(4)=1.5; x(5)=1.0; x(6)=0.7; stem(t,x); axis([-2.2 3.2 0 1.7]) 3、 单位阶跃函数 1,0()0,0t u t t ?≥?=?
常用函数 1 图形化信号处理工具，fdatool（滤波器设计），fvtool（图形化滤波器参数查看）sptool （信号处理），fvtool(b,a)，wintool窗函数设计.或者使用工具箱 filter design设计。当使用离散的福利叶变换方法分析频域中的信号时，傅里叶变换时可能引起漏谱，因此 需要采用平滑窗， 2数字滤波器和采样频率的关系。 如果一个数字滤波器的采样率为 FS，那么这个滤波器的分析带宽为Fs/2。也就是说这 个滤波器只可以分析[0，Fs/2]的信号.举个例字： 有两个信号，S1频率为20KHz，S2频率为40KHz，要通过数字方法滤除S2。 你的滤波器的采样率至少要为Fs=80HKz，否则就分析不到 S2了，更不可能将它滤掉 了！（当然根据采样定理，你的采样率 F0也必须大于80HK，,Fs和 F0之间没关系不大，可以任取，只要满足上述关系就行。） 3 两组数据的相关性分析 r=corrcoef(x,y) 4 expm 求矩阵的整体的 exp 4 离散快速傅里叶 fft信号处理中，傅里叶变换的典型用途是将信号分解成幅值分量和频率分量）。Ft为连续傅里叶变换。反傅里叶 ifft 5 ztrans（），Z变换是把离散的数字信号从时域转为频率 6 laplace（）拉普拉斯变换是把连续的的信号从时域转为频域 7 sound(x)会在音响里产生 x所对应的声音 8 norm求范数，det行列式，rank求秩 9 模拟频率，数字频率，模拟角频率关系 模拟频率f：每秒经历多少个周期，单位Hz，即1/s； 模拟角频率Ω是指每秒经历多少弧度，单位rad/s； 数字频率w：每个采样点间隔之间的弧度，单位rad。 Ω=2pi*f； w = Ω*T 10 RMS求法 Rms = sqrt(sum(P.^2))或者norm(x)/sqrt(length(x) var方差的开方是std标准差，RMS应该是norm(x)/sqrt(length(x))吧. 求矩阵的RMS：std(A(:)) 11 ftshift 作用：将零频点移到频谱的中间 12 filtfilt零相位滤波， 采用两次滤波消除系统的非线性相位， y = filtfilt(b,a,x)；注意x的长度必须是滤波器阶数的3倍以上，滤波器的 阶数由max(length(b)-1,length(a)-1)确定。 13 [h,t]=impz(b,a,n,fs)，计算滤波器的冲激响应 h为n点冲击响应向量 [h,x]=freqz(b,a,n,fs)计算频响，有fs时，x为频率f，无fs，x为w角频率， 常用于查看滤波器的频率特性 14 zplane(z,p) 画图零极点分布图 15 beta=unwarp(alpha) 相位会在穿越+-180发生回绕，可将回绕的 16 stepz 求数字滤波器的阶跃响应 [h,t] = stepz(b,a,n,fs) fvtool(b1,a1,b2,a2,...bn,an) fvtool(Hd1,Hd2,...) h = fvtool(...) 15 IIR数字滤波器设计方法 1 先根据已知带同参数求出最佳滤波器阶数和截止频率 [n,Wn] = buttord(Wp,Ws,Rp,Rs);
《数字信号处理》课程设计设计题目：基于MATLAB的音乐信号处理和分析 院系：物理工程学院 专业：电子信息科学与技术 学号： 姓名：
一、课程设计的目的 本课程设计通过对音乐信号的采样、抽取、调制解调、滤波、去噪等多种处理过程的理论分析和MATLAB实现，使学生进一步巩固数字信号处理的基本概念、理论以及频谱分析方法和数字滤波器设计方法；使学生掌握的基本理论和分析方法只是得到进一步扩展；使学生能有效地将理论和实际紧密结合；增强学生软件编程实现能力和解决实际问题的能力。 二、课程设计的基本要求 1 学会MATLAB的使用，掌握MATLAB的基本编程语句。 2 掌握在Windows环境下音乐信号采集的方法。 3 掌握数字信号处理的基本概念、基本理论和基本方法。 4 掌握MATLAB设计FIR和IIR数字滤波器的方法。 5 掌握使用MATLAB处理数字信号、进行频谱分析、涉及数字滤波器的编程方法。 三、课程设计内容 实验1音乐信号的音谱和频谱观察 使用windows下的录音机录制一段音乐信号或采用其它软件截取一段音乐信号（要求：时间不超过5s、文件格式为wav文件） ①使用wavread语句读取音乐信号，获取抽样率；（注意：读取的信号时双声道信号，即为双列向量，需要分列处理）； ②输出音乐信号的波形和频谱，观察现象； 使用sound语句播放音乐信号，注意不同抽样率下的音调变化，解释现象。 程序如下： [Y,FS,NBITS]=WAVREAD('怒放的生命 - 汪峰5s'); %读取音乐信号 plot(Y); %显示音乐信号的波形和频谱 sound(Y,FS); %听音乐（按照原来的抽样率） Y1=Y(:,1); %由双声道信号变为单声道信号 size(Y1) figure subplot(2,1,1);
南京航空航天大学基于Matlab的脑电信号处理 姓名陆想想 专业领域生物医学工程 课程名称数字信号处理 二О一三年四月
摘要：脑电信号属于非平稳随机信号，且易受到各种噪声干扰。本文基于Matlab仿真系统，主要研究了小波变换在脑电信号处理方面的应用，包括小波变换自动阈值去噪处理、强制去噪处理，以α波为例，提取小波分解得到的各层频率段的信号，并做了一定的分析和评价。关键词：脑电信号；小波变换；去噪重构；频谱分析 0引言 脑电信号EEG(Electroencephalograph)是人体一种基本生理信号，蕴涵着丰富的生理、心理及病理信息，脑电信号的分析及处理无论是在临床上对一些脑疾病的诊断和治疗，还是在脑认知科学研究领域都是十分重要的。由于脑电信号的非平稳性且极易受到各种噪声干扰，特别是工频干扰。因此消除原始脑电数据中的噪声，更好地获取反映大脑活动和状态的有用信息是进行脑电分析的一个重要前提。本文的研究目的是利用脑电采集仪器获得的脑电信号，利用Fourier变换、小波变换等方法对脑电信号进行分析处理，以提取脑电信号α波的“梭形”节律，并对脑电信号进行功率谱分析和去噪重构。 1实验原理和方法 1.1实验原理 1.1.1脑电信号 根据频率和振幅的不同，可以将脑电波分为4种基本类型[1]，即δ波、θ波、α波、β波。4种波形的起源和功能也不相同，如图1所示。 图1脑电图的四种基本波形 α波的频率为8~13Hz，振幅为为20~100μV，它是节律性脑电波中最明显的波，整个皮层均可产生α波。正常成人在清醒、安静、闭目时，波幅呈现有小变大，再由大变小，如此反复进行，形成所谓α节律的“梭形”。每一“梭形”持续时间约为1~2s。当被试者睁眼、警觉、思考问题或接受其他刺激时，α波立即消失而代之以快波，这种现象称之为“α波阻断”。一般
《数字信号处理》 （一）实验目的 使用stem绘图函数分别画出离散时间信号在指定范围内的图形。画图时使用xlabel,ylabel,title,legend等函数进行注释。复习MATLAB的基本应用，如：函数的定义、画图……并巩固理论知识中的多种离散时间信号及其图形，以及延迟与翻褶的函数变换等。（二）程序的运行与截图 1)用stem绘制单位阶跃序列u（n） clear all;close all;clc;%清除所有变量 n=0:50;%取值范围 y=(n>=0);%n>=0,y=1;n<0,y=0 stem(n,y);%显示出当0<=n<=50 时，函数u(n)的取值范围 xlabel('n');%对横轴进行注释 ylabel('y=u(n)');%对纵轴进行注释 title('y=u(n)的图形');%对图像的标题进行注释 legend('y=u(n)',2);%对图中曲线进行注释，标注在第二象限
2)用stem绘制单位抽样（冲激）序列δ（n）clear all;close all;clc; %清除所有变量 n=0:50; %取值范围 y=(n==0);%n=0,y=1;n!=0,y=1 stem(n,y);%显示出当0<=n<=50 时，函数δ(n)的取值范围 xlabel('n');%对横轴进行注释 ylabel('y=δ(n)');%对纵轴进行注释 title('y=δ(n)的图形');%对图像的标题进行注释 legend('y=δ(n)',2);%对图中曲线进行注释，标注在第二象限
3）用stem绘制矩形序列Rn（n）clear all;close all;clc; %清除所有变量 n=0:50; %取值范围 R10=((n>=0)&(n-9)<=0);%0<=n<=10,y=1;n>10,y=0 stem(n,R10);%显示出当0<=n<=50 时，函数Rn（n）的取值范围xlabel('n');%对横轴进行注释 ylabel(' y=R10（n）');%对纵轴进行注释 title('y=R10（n）的图形');%对图像的标题进行注释 legend('y=R10（n）',2);%对图中曲线进行注释，标注在第二象限