生物医学信号处理习题集
第一章 生物医学信号处理绪论 ..................................................................................................... 1 第二章 数字信号处理基础 ............................................................................................................. 1 第三章 随机信号基础 ..................................................................................................................... 5 第四章 数字卷积和数字相关 ......................................................................................................... 9 第五章 维纳滤波 ........................................................................................................................... 10 第六章 卡尔曼滤波 ....................................................................................................................... 13 第七章 参数模型 ........................................................................................................................... 16 第八章
自适应信号处理 (17)
第一章 生物医学信号处理绪论
1. 生物医学信号处理的对象是什么信号? 解答:
包括生理过程自发产生的信号,如心电、脑电、肌电、眼电、胃电等电生理信号和血压、体温、脉搏、呼吸等非电生理信号;还有外界施加于人体的被动信号,如超声波、同位素、X 射线等。
2. 生物信号的主要特点是什么? 解答:
随机性强,噪声背景强。
第二章 数字信号处理基础
You can use Matlab where you think it ’s appropriate. 1.FIR 滤波器和IIR 滤波器的主要区别是什么? 解答:
FIR 滤波器的单位脉冲响应是有限长的序列,该滤波器没有极点,具有稳定性。 IIR 滤波器的单位脉冲响应是无限长的序列,该滤波器有极点,有可能不稳定。
2.两个滤波器级联,第一个的传递函数为2-11z 2z 1)z (H -++=,第二个为-12z 1)z (H -=,当输入为单位脉冲时,求输出序列,画出级联滤波器的频率响应。 解答:
)z 1)(z 2z 1()z (H 12-1---++==32-1z z z 1----+
h(n)=[1,1,-1,-1],n=0,1,2,3。即输入单位脉冲时的输出序列值。 freqz(h,1)
3.A 3rd -order lowpass filter is described by the difference equation: )
3n (2781y .0)2n (1829y .1)1n (76y .1)
3n (0181x .0)2n (0543x .0)1n (0543x .0)n (0181x .0)n (y -+---+-+-+-+=
Plot the magnitude and the phase response of this filter and verify that it is a lowpass filter. 解答:
b = [0.0181, 0.0543, 0.0543, 0.0181]; a = [1.0000, -1.7600, 1.1829, -0.2781]; m = 0:length(b)-1; l = 0:length(a)-1; K = 500; k = 1:1:K;
w = pi*k/K; % [0, pi] 分成501个点. num = b * exp(-j*m'*w); % 分子计算 den = a * exp(-j*l'*w); % 分母计算 H = num ./ den;
magH = abs(H); angH = angle(H); subplot(1,1,1);
subplot(2,1,1); plot(w/pi,magH); grid; axis([0,1,0,1]) xlabel(''); ylabel('|H|'); title('幅度响应');
subplot(2,1,2); plot(w/pi,angH/pi); grid on; axis([0,1,-1,1])
xlabel('以pi 为单位的频率'); ylabel('以pi 弧度为单位的相位'); title('相位响应');
或freqz(b,a)
明显是低通滤波器,Wc 大概在0.25pi 。(衰减3个dB ,下降一半)
4.Find the inverse z-transform of x(z)=1
4z 3z z
2+-.To check the result using Matlab function residuez.
解答:
)z 3
111
z 11(2/1z 4z 3z 14z 3z z )z (X 1
12112--------=+-=+-=
b = [0,1]; a = [3,-4,1]; [R,p,C] = residuez(b,a) [b,a] = residuez(R,p,C) R = 0.5000 -0.5000 p = 1.0000 0.3333 C = []
b = -0.0000 0.3333
a = 1.0000 -1.3333 0.3333 笔算和程序结果一致。
5.Choose an appropriate window to design a digital FIR lowpass filter with the following specifications:
25dB .0A ,2.0p p ==πω,50dB A ,3.0s s ==πω
Determine the impulse response and provide a plot of the frequency response of the designed filter. (help fir1 function ) 解答:
wp = 0.2*pi; ws = 0.3*pi; tr_width = ws – wp ;
M = ceil(6.6*pi/tr_width) ;%查表求得窗长度,hamming window 即可 n=[0:1:M-1]; wc = (ws+wp)/2 b= fir1(M,wc/pi); h=b(1:end-1);
[hh,w] = freqz(h,[1],'whole');%默认就是hamming window hhh=hh(1:255);ww=w(1:255); % 画图
subplot(2,2,3); stem(n,h);title('实际脉冲响应') axis([0 M-1 -0.1 0.3]); xlabel('n'); ylabel('h(n)')
subplot(2,2,4); plot(ww/pi,20*log10(abs(hhh)));title('幅度响应(单位: dB )');grid axis([0 1 -100 10]); xlabel('频率(单位:pi )'); ylabel('分贝') set(gca,'XTickMode','manual','XTick',[0,0.2,0.3,1]) set(gca,'YTickMode','manual','YTick',[-50,0])
第三章 随机信号基础
1.什么是平稳各态遍历的随机过程? 解答:
如果随机信号的统计特性与开始进行统计分析的时刻无关,则为平稳随机过程,否则为非平稳随机过程。对于平稳过程,如果所有样本在固定时刻的统计特征和单一样本在全时间上的统计特征一致,则为各态遍历的随机过程。平稳且各态遍历是本课程分析医学信号的一个前提假设
2.判断随机相位正弦波在均值意义下是否各态遍历。)sin()(0φω+=t A t x ,A 0ω是固定值,φ是随机变量,分布为均匀分布:πφπ
φ20,21
)(<≤=p ,其它为零。 解答:
该随机过程的时间平均为:0dt )t (Asin 2T 1lim
m T
T
0T x =?+=-∞→φω 该随机过程的总体平均为:0d )(p )t (Asin dx )x (xp )x (E 20
0A
A
-=?+=?=φφφωπ 因此该过程在均值意义下是各态遍历的。
3.讨论相互独立、互不相关、相互正交的区别和联系。 解答:
随机变量统计独立的条件为:)y (p )x (p )y ,x (p = 互不相关的条件为:0)y ,x (cov = 正交的条件为:0)xy (E =
对于一般的随机变量:统计独立则互不相关;当其中有任意一个变量的均值为零,则互不相关和正交可以互相推导。
对于高斯随机变量,统计独立和互不相关可以相互推导;当其中有任意一个变量的均值为零,则三者都能互相推导。
4.输入序列n x 的一阶概率密度函数是
)
(2)(2n x n x u e x p n -=。证明:
5
.0)(=n x E ;如2142x x y +=,
1x 、2x 都是具有上述分布的随机序列,求)(y E 。
解答:
?=+∞
∞-n n n n dx )x (p x )x (E ?=+∞
0n -2x n dx 2e x n ?=+∞
0-2x n n de x -
?+∞+-=+∞-0
n 2x -2x n dx e 0e x n n
?=+∞0n 2x -dx e n ?-=+∞02x -n de 5.0=0)e (5.0n 2x ∞
+--=0.5 E (y )=E(2x 1+4x 2)=E(2x 1)+E(4x 2)=3
5.已知平稳随机过程x 的自相关函数如下,求其功率谱密度及均方,并根据所得结果说明该随机过程是否含有直流分量或周期性分量。
(ⅰ)πτπτττ
cos3cos 4e )(R x +=- (ⅱ)16cos 25e )(R 04x +=-τωττ
解答:
(ⅰ)ττωωτd e )(R )(P j x x -+∞
∞
-?=])(11
)(11[8)]3()3([22
πωπωπωδπωδπ-+++++-++= 514)0(R )x (E x 2=+==
因为0]12
[8)0(P 2
x ≠+=π,所以含有直流分量;
因为周期信号的自相关函数也是周期性的,而R 中包含有一个周期性的成分,因此该随机过程含有周期性分量。
(ⅱ)ττωωτd e )(R )(P j x x -+∞∞-?=])(161
)(161[50)(322020ωωωωωπδ-+++++=
411625)0(R )x (E x 2=+==
因为0]162
[
5032)0(P 2
x ≠++=ωπ,所以含有直流分量; 因为周期信号的自相关函数也是周期性的,而R 中没有包含周期性的成分,因此该随机过程不含有周期性分量。
6.设x (t )是平稳过程,)T t (x )t (x )t (y -+=,证明的)t (y 功率谱是:
)T cos 1)((2P )(P x y ωωω+=
解答:
ττωωτd e )(R )(P j y y -+∞
∞-?=Θ
其中))t (y )t (y (E )(R y ττ+=))]T t (x )t (x ())T t (x )t (x [(E -+++?-+=ττ )]T t (x )T t (x )T t (x )t (x )t (x )T t (x )t (x )t (x [E -+-+-+++-++=ττττ )(R )T (R )T (R )(R x x x x ττττ+-+++=
ττωωτd e )(R )(P j y y -+∞
∞-?=∴ττττωτd e )]T (R )T (R )([2R j x x x ?-+++=+∞
∞--
)(2P x ω=τττωτd e )]T (R )T ([R j x x ?-++++∞
∞-- ]e e 1)[(2P T j T j x ωωω++=-
)T cos 1)((2P x ωω+=,得证。
7.一个随机信号)t (x 1的自相关函数是τ
τ-=e A )(R 11,另一个随机信号)t (x 2的自相关函数为
τ
τ-=e
A )(R 22,在下列条件下,分别求信号相加后)t (x )t (x x(t)21+=的自相关函数)(R x τ。
(ⅰ))t (x 1,)t (x 2相互独立;
(ⅱ))t (x 1,)t (x 2来自同一信号源,只是幅度差一个常数因子K (K 不为1):)t (x 2=K )t (x 1。 解答:
(ⅰ))t (x 1,)t (x 2相互独立
))t (x )t (x (E )(R x ττ+=))]t (x )t (x ())t (x )t (x [(E 2121ττ+++?+= )]t (x )t (x )t (x )t (x )t (x )t (x )t (x )t (x [E 22211211ττττ+++++++= )(R )]t (x [E )]t (x [E )]t (x [E )]t (x [E )(R 221211ττττ+++++=
∞
→ττ)(limR 1Θ)]t (x [E ))t (x (E ))t (x (E )](t (t)x limE[x 121111=+?=+=∞
→τττ
0)]t (x [E 12=∴,同理0)]t (x [E 22=∴
)(R )(R )(R 21x τττ+=∴τ
-+=)e
A A (21
(ⅱ))t (x 2=K )t (x 1
由前面计算可得)(R x τ)(R )]t (x )t (x [E )]t (x )t (x [E )(R 221211ττττ+++++= )(R )]t (Kx )t (x [E )]t (Kx )t (x [E )(R 211111ττττ+++++= )(R )(2KR )(R 211τττ++= τ
-++=)e
A 2KA A (211
8.
n
x 是零均值,方差为2
x σ的白噪过程,把它先送入一个平均器,得)(2
1
1-+=
n n n x x y ,然后再将n y 送给一个差分器1y y z --=n n n ,求n z 的均值、方差、自相关函数和功率谱密度。 解答: 1y y z --=n n n )x x (21)(212n 1-n 1--+-+=
n n x x )(2
1
2--=n n x x 0))x x (2
1
(E )z (E 2n n n =-=-
)z (E 2n 2z =σ)z (E n 2-])(41[E 22--=n n x x )]x 2x (41[E 2
2n n 22-+-=-n n x x =2x σ/2=)0(R n z
)1(R n z =)z z (E 1n n +)]x x )((41
[E 1n 1n 2-+---=n n x x =0
)2(R n z =)z z (E 2n n +)]x x )((4
1[E n 2n 2--=+-n n x x =-2
x σ/4
当|m|>=3,自相关都为0。
)m (R n z =)]2m (21
)2m (21)m ([212x +---δδδσ
)w (P n z ]e 2
1e 211[21jw2jw2-2x --=
σ]cos2w 1[212x -=σ
9.随机序列
n
x 各次采样互相独立,且均匀分布于-1~1之间,设1--=n n n x x y ,
212--++=n n n n x x x z ,n n n x w w +-=-15.0,求n x 的均值和方差;n y 、n z 、n w 的自相关函数和
功率谱。
解答:
均值:?=+∞∞-n n n n dx )x (p x )x (E =?=+-1
1n n 00.5dx x
均方值:?=+∞∞-n n 22dx )x (p x )x (E n n 3
10.5dx x 11n 2
n
?==+- 方差:Var(x n )=3
1
)x (E )x (E n 2
2n =
- )0(R n y 2
12n )(E )y (E --==n n x x =
3
2
3
1
)]x x )(x x (E[)y y (E )1(R n 1n 1n n 1n n y n -=--==+-+
0)]x x )(x x (E[)y y (E )2(R 1n 2n 1n n 2n n y n =--==++-+,当|m|>=2时,y 的自相关函数都为零
)m (R n y =)]1m ()1m ()m (2[3
1+---δδδ,)w (P n y ]e 21e 211[32jw jw ---=]cosw 1[32
-=。
)0(R n z 2
2n 12n )x 2(E )z (E --++==n n x x =2
34
)]x 2x x )(x 2x x (E[)z z (E )1(R 1n n 1n 2n 1n n 1n n z n =
++++==-+--+ 3
1)]x 2x x )(x 2x x (E[)z z (E )2(R n 1n 2n 2n 1n n 2n n z n =++++==++--+, 当|m|>=3时,z 的自相关函数都为零
)m (R n z =)]2m ()2m ()1m (4)1m (4)m (6[31
++-+++-+δδδδδ,
)w (P n z ]e e 4e 4e 6[31j2w 2jw jw jw -++++=-cos2w 3
2
cosw 382++=。
)0(R n w 212n )5w .0(E )w (E --==n n x =)w (E 4
1
)w x (E )x (E 21n 1n n 2
n --+
-,由于中间一项为零,所以有 )0(R n w =
)x (E 3
4
2n =94, (0)5R .0]w )5w .0x (E[)w w (E )1(R n n w n n 1n 1n n w -=-==++
)0(R )5.0()1(5R .-0)]0.5w x (E[w )w w (E )2(R n n n w 2w 1n 2n n 2n n w -==-==+++,
所以)m (R n w =)0(R )5.0(n w m
-m
)5.0(9
4-=
, )(P n w ω)]m (R [DTFT n w =])5.0(94[DTFT m
-=ωω
ω
ωcos 4
53
1)5e .015e .05e .011(94j j j +=+-++=-。
第四章 数字卷积和数字相关
1.设x=[1 2 3 4];h=[4 3 2 1];求conv(x,h)、filter(h,1,x)、filter(x,1,h)的结果,并写出后两个函数对应的传递函数。 解答:
conv(x,h) =[4 11 20 30 20 11 4]
由于卷积的前后互换性,filter(h,1,x)=filter(x,1,h)=[4 11 20 30] H1(z )=321
234---+++z z z
,H2(z )=3214321---+++z z z
2.输入到线性系统的平稳随机过程x 是零均值、方差为1,输出信号的功率谱为
5.0,cos 211
)(2
=+-=
a a
a e S j y ωω,求此系统的传递函数H (Z )。要求该系统稳定。 解答:
,1)(=ωj x e S Θ
5.0,cos 211
)()()()(2
2
2=+-=
==a a
a e H e H e S e S j j j x j y ωωωωω, 因为系统要求稳定
1
5.011
)(,5.011)(---=
-=
∴z z H e e H j j ωω
3.列举相关技术在生物医学信号处理中的部分应用。 解答:
从噪声中检测信号,例如检测超声脉冲回波。
估计两个相似信号的时间延迟,例如测定微血管中的红血球流速,提取脑电诱发响应。 用于生物系统的辨识等。
4.估计相关函数时如果采用∑==+1
-m -N 0
n m
n n x x x N 1)m (R ?,估计是否为无偏估计?是否为一致估计? 解答:
是渐进无偏。
所以不是无偏估计,但,偏差了时,越大,偏差越大,偏差:
)N (R 0)]m (R ?E[N m m )m (R )]m (R
?E[,N )m (R N m N ]x E[x N 1)]m (R ?E[x x x
x
x 1
m N 0n m n n x ===∞→-∑=
=--=+
估计的方差当N 无穷时,趋于零。
因此该估计法是一致估计。
5.已知心电图的频率上限约为50Hz ,采集数据时候的采样频率至少为多少?如果采样频率为300Hz ,要求的频率分辨率为1Hz ,试确定做谱估计时每段数据的点数。 解答:
由于采样频率至少要为信号最高频率的两倍,所以这时采样频率至少要为100Hz 。
N f f s =
?Θ,300f
f
N s =?≥,做谱估计时每段数据的点数要大于300,考虑DFT 的计算,最好取2的幂次,可以取512点数据,或者补零到512点长。 第五章 维纳滤波
以下三题的系统模型图都参看该图。
1.设上图滤波器的方程是1-n n n x x y -=,输入s n 是确定性信号s n =b n ,b 是常数。n n 是白噪序列,零均
值,方差为2
n σ。求
(ⅰ)输出中的信号分量;
(ⅱ)输出中噪声分量的均方和方差; (ⅲ)输出中噪声分量的功率谱。 解答:
(ⅰ)输出中的信号分量
1n n 1n n 1n 1n n n 1-n n n n n b b n s n s x x y -----+-=--+=-=
前两项是由于信号s引起的输出,后面两项是噪声分量引起的输出。 (ⅱ)输出中噪声分量的均方和方差
0]n n [E 1n n =--,因此输出中噪声分量的均方等于方差:
2
n 1-n n 21-n 2n 21n n 2]n 2n n n [E ])n n [(E σ=-+=--
(ⅲ)输出中噪声分量的功率谱
由于1z 1)z (H --=,则ωωj j e 1)e (H --=,2
n j n )e (P σω=
)e (P )e (H )e (P j n 2
j j noise ωωω=∴)2cos 2(2
n
ωσ-=
2.对于上图中的系统模型,假设h (n )是因果的,用相关函数法推导出维纳滤波器的维纳-霍夫方程的离散形式,以及从该方程中解出了最佳滤波器)(n h opt 后的最小均方误差的最简式。 解答:
0,0)(<=n n h 当,∑+∞
=-==0)()()(?)(m m n x m h n s n y ,[
]
??
?
???--=∑+∞
=202
))()()(()(m m n x m h n s E n e E ,
Λ2,1,00
)())()()((20==??
????---∑+∞
=j j n x m n x m h n s E m opt
[][]
0)()()()()(0
≥--=-∑+∞
=j j n x m n x E m h j n x n s E m opt
维纳-霍夫方程:0)
()()(0
≥-=
∑+∞
=j m j R m h
j R m xx opt
xs
最小均方误差为:[
]
??
????--=∑+∞
=20min
2
))()()(()
(m opt m n x m h n s E n e E
])()()()()()()(2)([00
2
∑∑∑+∞=+∞
=+∞
=--+--=m r opt opt m r n x r h m n x m h m n x m h n s n s E
∑∑∑+∞
=+∞=+∞=??
?
???-+-=000)()()()()(2)0(m r xx opt opt m xs opt ss r m R r h m h m R m h R
∑+∞
=-=0min 2
)()()0()]([m xs opt ss m R m h R n e E
3.设线性系统如上图所示,已知n n n s ,相互独立,且ωω
2sin )(=j s e S ,2
1
)(=
ωj n e S 。要求设计一个滤波器ωω
2sin )(c e
H j =,试确定c 使得滤波后的输出n s
?与真实信号n s 的均方误差最小,即])?[(2n n s
s E -最小。 解答:
设误差为n n n s
?s e -= 其自相关为:
)m (R )m (R )m (R )m (R )]s ?s )(s ?s [(E )e e (E )m (R s ?s s ?s ?s s m n m n n n m n n e +--=--==+++
做傅立叶变化:)()()()()(???ωωωωω
j s j s s j s s j s j e e S e S e S e S e
S +--=
ωωωωωωωω4262j n j s 2j j x 2j ?sin 2
1
sin ])(e S )(e S [)e (H )(e S )e (H )(c c e S j s +=+==
ωωωωωω4i s i i sx i ?sin )e (S )e (H )e (S )e (H )(c e S j s s === ωωωωωω4i s i i xs i s ?sin )e (S )e (H )e (S )e (H )(c e S j s ===**
22
14321
c c +-=ξ
求导等于零:4
3=opt c
4.)n (x 是白噪过程,零均值,方差为2
x σ,把)n (x 作为输入加到一线性系统上,系统的冲激响应是
)n (h ,输出是)n (y 。证明:
(ⅰ))0(h )]n (y )n (x [E 2
x
σ= (ⅱ)∑=+∞
∞
=-n 2
2x 2y )n (h σσ
解答:
(ⅰ)∑=∑-=∑-=m
x m m
)
m (R )m (h )]
m n (x )n (x [E )m (h ]
)m (h )m n (x )n (x [E )]n (y )n (x [E
当m=0时,2
x x )0(R σ=,当m 为其它时0)m (R x =,代入上式,
)0(h )]n (y )n (x [E 2x σ=∴
(ⅱ)])k n (x )k (h )m n (x )m (h [E ]y y [E m
k
n n 2y ∑∑--==σ
∑∑--=m k
)]k n (x )m n (x [E )k (h )m (h
∑∑-=m k
x )k m (R )k (h )m (h
当m=k 时,2x x )k m (R σ=-,当m 为其它时0)k -m (R x =,代入上式,
∑=∑=+∞
∞
=-n 22x 2x k
22y )n (h (k)h σσσ
得证。
第六章 卡尔曼滤波
1.有一信号)n (s ,其自相关函数Λ2,1,0,7.0)(±±==m m R m
s ,被一零均值,方差为0.4的白噪)n (n 所淹没,)n (s 与)n (n 统计独立。
(ⅰ)设计一个长度等于3的FIR 数字滤波器,其输出)n (y 使得]))n (s y(n)
[(E 2-最小化。 (ⅱ)设计一个因果的最优滤波器,并说明如何在计算机上实现。
解答:
(ⅰ)根据均方误差最小准则得到W-H 方程:
1,,2,1,0)
()()(1
-=-=∑-=N j m j R m h j R N m xx opt xs Λ,其中x=s+n ,表示输入信号,
因为N=3,且)m (R )m (R )m (R n s xx +=,
)m (R )]m n (s )n(n)s(n)[(E )]m n (s )n (x [E )m (R s xs =++=+=,代入W-H 方程得: 2,1,0)]
()()[()(n 2
s s =-+∑-==j m j R m j R m h j R m opt
把Λ2,1,0,7.0)(±±==m m R m
s ,)m (4.0)m (R n δ=代入上式得三个方程:
???
???
?∑+==∑+==∑+=====2
0m m -2opt 220
m m
-1opt 20
m m -opt m)]
(20.4(m)[0.7h 0.7:2j m)](10.4(m)[0.7h 0.7:1j m)]
(0.4(m)[0.7h 1:0j ---δδδ ??
??
??????=??????????2opt opt opt 2
27.07.01)2(h )1(h )0(h 4.17.07.07.04.17.07.07.04.1 解得:?
????
?????=??????????0517.01681.00.6121)2(h )1(h )0(h opt opt opt 所以设计的滤波器的传递函数为:-210517z .01681z .06121.0)z (H ++=- (ⅱ)设计一个因果的最优滤波器 因为)m (R )m (R )m (R n s xx += 所以输入信号的z 变换为: 4.0)
7z .01)(7z .01(51
.0)z (R )z (R )z (R 1n s xx +--=
+=-
)
7z .01)(7z .01()az 1)(az 1(b )7z .01)(7z .01(0.28z -0.28z -1.10601
11-1----=--=---
列出方程求系数a 与b ,?
??==+28.0ab 1060
.128.a0b 利用solve 函数解出a ,b :
a=[3.68 ,0.2718],b=[0.76, 1.03],a 取小于1的数,所以a=0.2718,b=1.03
因此)7z .01)(7z .01()
2718z .01)(2718z .01(1.03)z (R 11xx ----=--
则03.12
w
=σ,)7z .01()2718z .01()z (B 1
1----=,)
7z .01()
2718z .01()z (B 1---= )(z H opt )()](/)([2
11z B z B z R w xs σ+-=)()](/)([211
z B z B z R w s σ+-=+---??????----=)2718.01)(7.01(51
.0)2718.01(03.17.01111z z z z ??????---=---)7.01(6298.0)2718.01(03.17.01111z z z ??????-=-)2718.01(03.16298
.01z )n (u )2718.0(6115.0)n (h n opt =
计算机实现可以利用均方误差∑-=-
=1
min 2
)()()0()]([N m ss opt
ss m R m h
R n e E ,当取N 与N+1时它们的
均方误差非常接近时就可以确定N 了。
2.比较维纳滤波和卡尔曼滤波方法的区别和联系。 解答:
维纳滤波和卡尔曼滤波都是解决线性滤波和预测问题的方法,并且都是以均方误差最小为准则的,在平稳条件下两者的稳态结果是一致的。但是它们解决问题的方法有很大区别。维纳滤波是根据全部过去观测值和当前观测值来估计信号的当前值,因此它的解形式是系统的传递函数)(z H 或单位脉冲响应)(n h ;卡尔曼滤波是用当前一个估计值和最近一个观测值来估计信号的当前值,它的解形式是状态变量值。维纳滤波只适用于平稳随机过程,卡尔曼滤波就没有这个限制。设计维纳滤波器要求已知信号与噪声的相关函数,设计卡尔曼滤波要求已知状态方程和量测方程,当然两者之间也有联系。
3.写出卡尔曼滤波的状态方程与量测方程,并解释。画出卡尔曼滤波的信号模型。 解答:
状态方程:1)(k w 1)A(k)S(k S(k)1-+-=
表示的含义就是在k 时刻的状态S(k)可以由它的前一个时刻的状态1)S(k -来求得,即认为k -1时刻以前的各状态都已记忆在状态1)S(k -中了。
量测方程:w(k)C(k)S(k)
X(k)+= 用X(k)表示量测到的信号矢量序列,w(k)表示量测时引入的误差矢量,C(k)称为量测矩阵,它的
引入原因是,量测矢量X(k)的维数不一定与状态矢量S(k)的维数相同,因为我们不一定能观测到所有需要的状态参数。 卡尔曼滤波的信号模型:
4.根据卡尔曼滤波的状态方程和量测方程,假设A(k)和C(k)是已知的,X(k)是观测到的数据,也
是已知的,假设信号的上一个估计值1)(k S
-?已知,如何来求当前时刻的估计值(k)S ?? 解答:
1)(k w 1)A(k)S(k S(k)1-+-=,w(k)C(k)S(k)X(k)+=
假设暂不考虑(k)w 1与w(k),用上两式得到的(k)S
?和(k)X ?分别用(k)S '?和(k)X '?表示,得: 1)(k S A(k)(k)S
-='??,1)(k S C(k)A(k)(k)S C(k)(k)X -='='??? 必然,观测值X(k)和估计值(k)X '?之间有误差,它们之间的差(k)X ~
称为新息(innovation ):
(k)X
X(k)(k)X '-=?~ 新息的产生是由于我们前面忽略了(k)w 1与w(k)所引起的,也就是说新息里面包含了(k)w 1与
w(k)的信息成分。因而我们用新息(k)X ~
乘以一个修正矩阵H(k),用它来代替式的(k)w 1来对S(k)
进行估计:
(k)X H(k)1)(k S A(k)(k)S ~??+-=
1)](k S C(k)A(k)H(K)[X(k)1)(k S
A(k)--+-=?? 1)](k S C(k)A(k)w (k)(k)H(K)[C(k)S 1)(k S A(k)(k)S
--++-=??? 1)](k S C(k)A(k)w (k)1)](k w 1)(k S A(k)H (K )[C(k)[1)(k S A(k)1
--+-+-+-=??? H (k)w (k)1)](k w 1)(k S (k)H (K )C(k)[A H (k)C(k)]1)[I (k S A(k)1
+-+-+--=?? 根据上式来求最小均方误差下的H(k),然后把求到的H(k)代回去则可以得到估计值(k)S
?。 ⊕
⊕
w(k)
(k)
w 1X(k)
第七章 参数模型
1.对于一个随机信号,可以对它进行频谱分析,叙述AR 谱法和周期图法相比的优点。 解答:
平滑、需要较短数据即可、频率分辨率高、峰值包络线的好估计等。
2.设已知Λ2,1,0),4.0(11
3)8.0(1114)(±±=-=
m m R m
m s ,用L-D 算法为此信号估计p =1,2,3阶AR 模型的系数和激励白噪的功率。
解答:
计算自相关函数:111
31114)0(=-=
s R 9091.04.0113
8.01114)1(=-=s R 7709.016.011
3
64.01114)2(=-=s R 6342.0064.011
3
512.01114)3(=-=
s R , 下面为了简写,省略下标s 。
按照L-D 算法得初始功率和系数为:1a 1,)0(E 00===s R P=1:
1736
.0)]1(a 1[E E 9091
.0E /)1(R )1(a 2
01011=-=-=-=
P=2:1558.0)]2(a 1[E E 2.1)]2(a 1)[1(a )1(a 32.0E )1()1(a )2()2(a 2212212
112=-=????
?
-=+==--=
R R
P=3:1558
.0)]2(a 1[E E 0.32(1)(3)a a (2)a )2(a 2.1)2((3)a a )1(a )1(a 0E )1(R )2(a )2()1(a )3()3(a 2212232323232223=-=???
????=+=-=+=≈---=R R 3.某随机过程用AR 模型拟和的结果是4
32158.06.258.45.311
)(----++++=
z z z z z H ,试由它导
出一个ARMA (2,1)模型。 解答:
设22111
1ARMA z a z a 1z b 1)z (A )z (B )z (H ---+++==,要使得它与H (z )相等则有 221111z a z a 1z b 1---+++4
32158.06.258.45.311
----++++=z
z z z
列出各系数方程: 0
58b .006b .258.0058b .46.2a 58.45.3b a 5.3b 1112
111==+=+=+=+,利用后三个方程,最小二乘法解出b 1来,然后再利用前两个方程求出a 1,a 2 ????
??????=??????????+??????????3211e e e 058.06.2b 58.06.258.4,两边同乘以[4.58 2.6 0.58] 得伪逆解b 1=-0.4779,所以a 1=3.0221,a 2=2.9074
ARMA (2,1)模型为:2
11ARMA 9074z
.20221z .310.4779z 1)z (H ---++-= 第八章 自适应信号处理
1.画出自适应噪声抵消的框图,并证明滤波后的输出将在最小均方意义下抵消噪声,同时,抵消后的结果将在最小均方意义下逼进信号。 解答:
因为2
2
))()(())()((n n n y E n e n s E -=-
))()((2))(())(())()((222n e n s E n e E n s E n e n s E -+=- ))()((2))()((2))(())((22n s n y E n n n s E n e E n s E +-+-= 22))(())((n e E n s E +-=
所以当均方误差最小时,s 与e 最逼近,同时y 与n 也最逼近。即滤波后的输出y 将在最小均方意义下抵消噪声n ,同时,抵消后的结果e 将在最小均方意义下逼进信号s 。
2.列举自适应信号处理在生物医学信号处理中的部分应用。 解答:
自适应噪声抵消:例如母腹电极上胎儿心电的提取,心电图中工频干扰的抑制,心电图中高频手术刀干扰的去除,呼吸阻抗中心电伪迹的消除等。 自适应谱线增强:诱发脑电的提取等。
自适应系统辨识:血压与血管收缩药物之间的系统自适应控制等。
3.简述横向结构的随机梯度法算法步骤。 解答:
步骤1:]x ,x ,x ,x [)T (X p 1p -T 2-T 1-T T '=+Λ?
个值观察到
步骤2:先给定。预先给出,与,初值计算μμT T e W (T)X e 2(T)W 1)(T W ?
??+=+ 步骤3:1)
(T X )1(T W d e ]x ,x ,x [)1T (X x 1T 1T 2p T T 1T 1T +'+-='
=++++-++?
?Λ?
计算新的误差:后,令当有新观测值 转入步骤2,代入得到W (T+2),e (T+2)…..使得W 接近最优解。停止循环的判断规则多样。