现代数字信号处理-使用MATLAB分析与实现（新形态版）习题及答案汇 第1-8章

1-1已知系统H(z)的单位采样响应hn=δn解：H(z)=1-因果离散时间LTI系统存在稳定因果逆系统的充分必要条件是：H(z)的所有零点（极点）必须在单位圆内（最小相位系统）零点：1-z=两个零点都在单位圆内，因此系统为最小相位系统，存在因果稳定的逆系统。逆系统G(z)=G(z)=其单位采样响应为：g(n)=或写为：g(n)=1-2一个线性移不变系统的传输函数为Hz=z（1）试求实现这个系统的差分方程。（2）试证明这个系统是一个全通系统。（3）H(z)和另一个系统G(z)级联后，整个系统函数为1，如果G(z)是一个稳定系统，试求其单位采样响应解：（1）H(z)=

Y(z)(1-a

y(n)-ay(n-1)=x(n-1)-因此差分方程为：

y(n)=ay(n-1)+x(n-1)-（2）

H(分子项：

e并且对任意复数

β

有

1-β取β=a*因此频率响应分子等于分母，等于1，得证。（3）H(z)G(z)=1

G(z)=

g(n)=1-3证明如下三个系统具有相同的幅频响应，并判断哪一个是最小相位系统，哪一个是最大相位系统，哪一个是混合相位系统。解:

HHH三个系统的幅频响应相同，因为它们的零极点模乘积比例相同最小相位系统：所有零点在单位圆内最大相位系统：所有零点在单位圆外混合相位系统：零点在单位圆内外都有1-4假设Hz和G（1）HzG(z)解：（1）已知Hz和G零极点分布：HzG(z)的零点是H因为H和G的所有零极点都在单位圆内，所以HzG因果稳定性：若H和G均因果稳定，则Hz因此Hz两个最小相位系统的和不一定是最小相位，因为零点会发生变化，可能出现在单位圆外。1-5一个因果线性移不变系统的传输函数为Hz=(1-2z-2)(解：

HHz1-6试用初值定理证明，如果hmn是最小相位序列，且hn解：设hmn是最小相位因果序列，则它们的z变换之间只相差一个全通函数HapH(z)=由初值定理（对因果序列）：h设最小相位系统函数为H当z∞时，lim取K>0（可通过相位适当选择使hm全通函数为H当zlim于是lim因为pk<1（若于是h(0)=取绝对值：h(0)由于c<1且hh(0)1-7表1.7.1所列的所有系统都是同态系统，已知其输入运算，试确定其输出运算。表1.7.1题1-7表序号系统变换T[x(n)]输入运算1y加法2y乘法3x加法4x卷积5x乘法6y乘法7y乘法8y加法9y乘法解：（1）输入加法：x输出：2所以输出运输也是加法。（2）T而T不相等，所以不是输入乘法。（3）x(n)=Z输出运算：加法（在频域上，两个z变换的加法）。（4）Z所以输出运算是乘法。（5）时域乘法zZ所以输出运算：复卷积（循环卷积）。（6）检查：

TT两者相等，所以输出运输也是乘法。（7）取绝对值：TT相等，所以输出运算是乘法。（8）TT所以输出运算：乘法。（9）TT与ex通常表述为：输出运算是对数域的乘法映射回原始域的运算，简称为幂运算(⨂)使得a⨂b=exp(lna∙lnb)。序号系统变换T[x(n)]输入运算输出运算1y加法加法2y乘法与第一行相同，应为加法3x加法加法4x卷积乘法5x乘法复卷积（或圆周卷积）6y乘法乘法7y乘法乘法8y加法乘法9y乘法对数域乘法（即⨂：a⨂b=exp(lna∙lnb)）1-8考虑通信系统中的信道均衡场景，接收信号y(n)=h(n)解：（1）数学模型基于二阶统计量的盲反卷积，接收信号:y其中(n)为独立信号，h(n)为未知FIR信道，v（2）LMS自适应反卷积滤波器设计滤波器输出：s其中ω(n)为均衡器系数，Y(n)=[y(n),y(n-1),...,y(n-M)]恒模准则（CMA）：代价函数：J(瞬时梯度下降（LMS型更新）：ω(n+1)=（3）迭代公式滤波：s误差：e(n)=更新：ω(n+1)=（4）步长

μ

的影响收敛速度：时间常数τ∝1(μ稳定性：需满足0<μ<2稳态失调：失调量Mdis1-9已知卷积同态系统的输入x(n)解：卷积同态系统的一般结构为：y(n)=其中D*L为线性系统(对复倒谱序列进行处理)。D*-1当x(n)=δX(z)=1复倒谱x(n)=0(对所有n于是：D若L是线性系统，那么L[0]=0。因此:y(n)=而逆特征系统D*-1对零序列的作用：频域Y(z)=0，取指数得y(n)=1-10什么是倒谱？什么是复倒谱？复倒谱有哪些主要性质？（1）倒谱对信号x(n)的幅度谱取对数，再求傅里叶逆变换得到的序列，称为实倒谱(通常简称为倒谱)定义(实倒谱)：c特点：只保留幅度信息，丢失相位信息。（2）复倒谱对信号x(n)的z变换取复对数，再进行逆z变换(或傅里叶逆变换)得到的序列，保留幅度与相位信息。定义：x其中logX(（3）复倒谱有哪些主要性质卷积性质：若x(n)=h(n)*s(n)，则x(n)=能量集中性：最小相位信号的复倒谱集中在n≥0，最大相位信号的复倒谱集中在指数/冲击序列关系：若x(n)稳定性：复倒谱x(n)当|n|对称性：对实信号，复倒谐满足x(n)=x递推计算：可通过递推公式从x(n)直接计算复倒谱，适用于最小相位信号(无需求解复对数相位展开2-1设一个随机信号为，其中为常数，A是随机变量，符合均值方差为的高斯分布，求均值和自相关函数，判断该信号是否为宽平稳的。解：均值：E[x(n)]=0自相关函数：r平稳性判断：该信号不是宽平稳信号，因为自相关函数与2-2对于平稳随机信号和，如果，证明解：由定义推导：r因此：r2-3若随机信号向量满足如下的高斯分布： 证明： 证明期望：E证明协方差矩阵：E2-4设离散时间过程按下式产生：，其中是方差为的白噪声过程。另一个过程是和噪声之和，即，是方差为的白噪声，且与不相关。试求：（1）的功率谱；（2）的功率谱。解：系统函数：H(1)

x(n)

的功率谱：P(2)

z(n)

的功率谱：P2-5设给定一个线性移不变系统，其系统函数为，它受零均值的指数相关噪声的激励产生随机过程。已知的自相关序列为，试求（1）的功率谱；（2）的自相关序列；（3）和之间的互相关；（4）互功率谱，它是互相关的Z变换。解：（1）P首先求

Px​(z)：P因此：P化简得：P（2）rr（3）r具体计算：r最终r（4）P代入得：P2-6考虑随机过程，其中是均值为零、方差为的高斯白噪声，试针对如下情况，求的自相关序列和功率谱：（1）A是均值为0，方差为的高斯随机变量，和是常量；（2）是在区间上均匀分布的随机变量，A和是常量；（3）是在区间上均匀分布的随机变量，A和是常量。解：已知：x（1）自相关序列：rEr该信号非平稳，自相关函数依赖于时间（2）此时信号平稳。自相关序列：r功率谱：P（3）自相关序列：r=E因此：r功率谱：P2-7确定如下的自相关阵是否合法？若不合法，请说明原因。（1）（2）（3）（4）（5）解：(1)R(2)R(3)R(4)R(5)R原因：对角元2-8随机过程，其中、为常数，为均匀分布于之间的随机变量，即 （1）求随机过程的自相关函数；（2）判断是否是平稳过程。解：已知：X自相关函数R利用三角恒等式：cosα-β=ωR第一项与

Θ

无关：a第二项：E===因此：R令

τ=t2​−t1​：R(2)平稳性判断均值：=均值依赖于

t，不是常数。因此：不平稳（均值随时间变化）2-9已知非平稳随机过程的自相关函数，求其功率谱密度。解：由于自相关函数只依赖于t，不依赖于t′，该过程为非平稳过程。对于非平稳过程，功率谱密度定义为自相关函数关于τ=t′−t

的傅里叶变换（时变功率谱）：S代入

R(t,t+τ)=21​(t+cosω0​t)（与

τ

无关）：S-因此：S2-10某随机过程由下述三个样本函数组成，且等概率发生： （1）计算数学期望和自相关函数；（2）该随机过程是否平稳？解：已知：X数学期望自相关函数R=因为sint1​sint2​+cost1​cost2​=cos(t1​−t2​)。令τ=t2​−t1​，则t1​−t2​=−τ，cos(−τ)=cosτ：R平稳性判断mXt=3-1设随机过程xn=i=1Kaicosωin+φi+w(n)，其中，a解：对于单个余弦项，其自相关为E[由于φiR对RxS3-2已知x(n)满足AR(2)模型xn+a1xn-1+a2解：自相关函数满足Yule-Walker方程：R根据观测数据x(0),…,x(N-1)估计出Rx(0),Rxa3-3二阶AR(2)过程xnx其中，w(n)为均值为0方差为0.5的白噪声。（1）写出该随机过程的Yule-Walker方程；（2）求xn解：（1）Yule-Walker方程根据AR(2)模型的差分方程，可知a1=0,aRRR(2)+0⋅R(1)+0.78R(0)=0解得R(0)=RR（2）Var3-4设平稳随机过程xn是由零均值方差为1的白噪声w(n)激励系统函数为HS试求系统函数Hz解：设H(z)最小相位，σwS=有一极点处于单位圆外。为了获得最小相位，需把单位圆外的极点用其倒数代替，这样S设z=eSxS取全部极零点在单位圆内H(z)，则H(z)=3-5已知平稳随机信号xn的自相关函数值R0=1,R1=0.5,R2=0.5,R3=0.25。现用3阶全极点模型（AR(3)模型）估计它的功率谱，设模型参数b0解：初始化，σ0p=1：aσp=2：aaσp=3：aaa由于b(0)≠1，则AR(3)系统差分方程为k=03k=0当m=0，R(0)+代入可得b(0)2σw2=213-6已知某自回归过程的5个观测值为{1,-1,0,-1,1}。（1）利用L-D算法设计一个3阶AR模型，确定模型的参数；（2）利用Burg算法求一阶和二阶反射系数，并画出二阶预测误差滤波器的格型结构图；（3）求该自回归过程的功率谱估计。解：（1）自相关函数R=利用L-D算法可以递推求出3阶AR模型的参数为aAR(3)模型为x(n)=-0.5（2）初始条件：f0n=b0k1阶前、后项预测误差：fb当p=2时k二阶预测误差滤波器的格型结构图为（3）AR(3)模型的功率谱：S(ω)=代入L-D结果，得S(ω)=3-7已知序列xn={4.684,7.247,8.428,8.650,8.640,8.392,7.321}由模型xn=1.7xn-1-0.72xn-2解：xn有非零均值，均值μ≈(4.684+7.247+8.428+8.650+8.640+8.392+7.321)/7≈7.623，对xy初始化fσ当p=1时kσ前后向误差fb代入，可得fb当p=2时kaa2σ与真实σw注意，这里估计的a1,a3-8一个自回归过程的5个观测值为x（1）用L-D算法设计一个二阶线性预测器，计算各阶线性预测系数、预测误差功率，以及x0和x4的预测值x0（2）用Burg算法计算一阶和二阶反射系数、预测误差功率，以及预测值x0和x4（3）对用上面两种算法求得的预测误差功率、预测值进行比较，看看哪种算法预测更准确些。解：（1）AR模型假设零均值过程，否则Yule-Walker方程不成立。因此，首先对xn去均值，得到去均值后yR可得RRR利用L-D算法，当p=1时σ当p=2时ka对y(n)，二阶y对于y(0)，需要y(-1),y(-2)，假设为0，所以y(0)=0，对于y(4)，预测值y4=-横向结构预测误差滤波器e（2）Burg算法对去均值后数据y(n)，初始f0kkσσyxx二阶格型结构预测误差滤波器的框图（3）在预测精度上，Burg预测x(4)更接近真实值5，误差更小。在预测误差功率上，Burg的σ22远小于L-D，说明模型拟合更好。3-9随机过程xnx其中，a和ω0是常数且a>0，φ是相互独立并在[0,2π]服从均匀分布的随机相位，w(n)为均值为0方差为σw2R试用Pisarenko谐波分解法确定正弦波的频率ω0、幅度a和白噪声的方差σ解：已知Rx(0)=3,RxR因此RRR求解可得ω3-10设随机信号x(n)是由加性白噪声和两个复正弦信号组成，其三阶自相关矩阵R（1）确定噪声的平均功率；（2）求信号子空间和噪声子空间；（3）判断信号子空间和噪声子空间是否正交；（4）用Pisarenko谐波分解法估计信号频率；（5）用MUSIC方法估计信号频率。解：（1）Rσ（2）信号子空间ES=span{v1噪声子空间EN=span{（3）因为v1Tv3=（4）f(z)=令f(z)=0，可得z=±j，因此ω1=（5）P令-12-14-1维纳滤波器的输出是否一定是“更清晰”的信号？为什么？‌‌解析：‌不一定‌。维纳滤波器输出的是‌在均方误差意义下最优的线性估计‌，但并不保证“视觉上更清晰”或“细节更丰富”。‌原因包括‌：若噪声估计偏高，会过度平滑，丢失细节（如边缘模糊）；若信号功率谱估计不准，可能增强噪声而非抑制；它是线性滤波，无法恢复非线性失真或超分辨信息。因此，维纳滤波器是“统计最优”，不一定是“感知最优”。‌4-2请说明维纳滤波器与低通滤波器在去噪应用中的本质区别。‌‌解析：‌对比项维纳滤波器低通滤波器设计依据基于信号与噪声的‌功率谱统计特性‌基于‌人为设定的截止频率‌自适应性自适应频段增益（信号强则通过，噪声强则抑制）固定增益，所有高频一律衰减信息利用利用信号与噪声的先验统计信息仅利用频率分布的直观假设效果在已知统计特性时更优，保留更多有效信号简单但易损失有用高频成分‌4-3为什么维纳滤波器被称为“最优线性滤波器”？这里的“最优”指什么？‌‌解析：‌维纳滤波器被称为“最优线性滤波器”，是因为在所有‌线性时不变系统‌中，它能‌最小化输出信号与原始理想信号之间的均方误差（MSE）‌。这里的“最优”是指：（1）在限定类别内最优‌：仅限于线性系统；‌（2）在统计意义下最优‌：基于最小均方误差准则；‌（3）有解析解‌：可通过傅里叶变换直接求得，无需迭代。它不是全局最优（非线性滤波器可能更好），但在线性框架下是理论最优解。‌4-4卡尔曼滤波器的核心思想是什么？‌‌解析‌：卡尔曼滤波器的核心思想是‌融合预测与观测‌，通过系统模型预测状态，并利用观测数据修正预测，从而在不确定性中找到最优估计。它基于线性高斯系统假设，通过预测步和更新步交替进行，实现状态估计的最优性。‌4-5卡尔曼增益K的作用是什么？‌‌解析‌：卡尔曼增益K的作用是‌决定预测值和测量值在更新中的权重‌。4-6已知观测信号，其中为为AR(1)过程：，为白噪声，方差，为白噪声，方差。求一阶维纳滤波器系数。解析：AR(1)过程自相关：观测信号自相关：维纳-霍夫方程：即：解方程组（用矩阵求逆或代入）：第一式乘1.36，第二式乘0.6：→→相减（第一式减第二式）：→代入第一式：→→最终结果：4-7已知离散时间系统中，期望信号d(n)与输入信号x(n)满足以下关系：输入信号x(n)=s(n)+v(n)，其中s(n)为有用信号，v(n)为加性白噪声；有用信号s(n)的自相关函数rss噪声v(n)与s(n采用长度为N=2的FIR维纳滤波器，即滤波器输出y(n)=w0试求：输入信号x(n)的自相关函数rxx期望信号d(n)与输入信号x(n)的互相关函数rdx最优滤波器系数w0*和解析：维纳滤波器的最优系数满足维纳-霍夫方程：Rxx其中：Rxxw*rStep1：计算输入信号的自相关函数r由于x(n)=s(n)+v(n)，且v(n)分别计算k=0和k=1的值：当k=0时：rss(0)=0.8故r当k=1时：rss(1)=故rxx（注：自相关函数为偶函数，rxxStep2：计算互相关函数r题目中期望信号d(n)本质是“恢复有用信号s(n)”，因此d(n)=s(n)（默认维纳滤波的目标是最小化E[互相关函数定义为rdx(kr由于v(n)与srdx分别计算k=0和k=1的值：rdx0=rStep3：构建维纳-霍夫方程并求解最优系数对于N=2的FIR滤波器，自相关矩阵RxxR互相关向量维纳-霍夫方程为：1.2将方程展开为线性方程组：1.2w0.8w求解该方程组：式1×1.2：1.44w式2×0.8：0.64w式3-式4：0.8将w0*=0.7最终结果:输入自相关函数：rxx互相关函数：rdx最优滤波器系数：w4-8假设一个目标沿直线匀速运动，仅通过传感器测量其位置（无速度直接观测）。系统模型如下：状态向量：其中为位置，为速度。状态转移矩阵（采样周期秒）：观测矩阵：过程噪声协方差：观测噪声方差：初始状态估计：实际观测序列（位置测量值）：请计算第1步（k=1）的预测与更新结果，即和。解析：Step1:预测(k=1)先计算中间乘积：再乘：加Q:Step2:更新(k=1)计算卡尔曼增益：先算:加R：计算状态更新：计算协方差更新：最终答案(k=1):4-9假设无人机使用GPS（测位置）和IMU（测速度）进行融合定位。状态向量：系统模型（采样周期T=0.1s）：观测模型（同时观测位置和速度）：初始估计：真实观测值(k=1)：请计算第1步的卡尔曼增益和更新后的状态估计。解析：Step1:预测先算：再乘：加Q：Step2:更新由于H=I，卡尔曼增益简化为：先算+:求逆（矩阵公式）：设行列式：逆矩阵：]计算逐项计算：第一行第一列：第一行第二列：第二行第一列：第二行第二列：所以：更新状态：最终答案：4-10通过MATLAB的wiener2函数，实现维纳滤波器对添加了高斯噪声的经典图像的复原。解析：MATLAB代码示例‌：%读取图像并添加高斯噪声original_img=imread('cameraman.tif');noisy_img=imnoise(original_img,'gaussian',0,0.01);%使用wiener2进行维纳滤波（局部自适应）filtered_img=wiener2(noisy_img,[55]);%[55]为邻域窗口大小%显示结果对比figure;subplot(1,3,1),imshow(original_img),title('原图');subplot(1,3,2),imshow(noisy_img),title('加噪图像');subplot(1,3,3),imshow(filtered_img),title('维纳滤波结果');演示结果如下：6.MATLAB中实现卡尔曼滤波的代码通常包括状态预测、协方差预测、卡尔曼增益计算、状态更新和协方差更新等步骤。试写一个简单的MATLAB代码示例，演示如何实现卡尔曼滤波。解析：MATLAB代码示例：%初始化参数A=1;%状态转移矩阵B=0;%控制输入矩阵H=1;%测量矩阵Q=0.1;%过程噪声协方差R=1;%测量噪声协方差x0=0;%初始状态P0=1;%初始误差协方差dt=1;%时间步长T=20;%总时间time=0:dt:T;%时间向量u=zeros(size(time));%控制输入%真实状态x_true=zeros(size(time));x_true(1)=x0;fort=2:length(time)x_true(t)=A*x_true(t-1)+sqrt(Q)*randn;end%测量数据y_meas=H*x_true+sqrt(R)*randn(size(time));%卡尔曼滤波估计x_est=zeros(size(time));x_est(1)=x0;P=P0;fort=2:length(time)%预测x_pred=A*x_est(t-1)+B*u(t);P_pred=A*P*A'+Q;%更新K=P_pred*H'/(H*P_pred*H'+R);x_est(t)=x_pred+K*(y_meas(t)-H*x_pred);P=(1-K*H)*P_pred;end%绘制结果figure;plot(time,x_true,'g','DisplayName','TrueState');holdon;plot(time,y_meas,'r.','DisplayName','Measurements');plot(time,x_est,'b--','DisplayName','EstimatedState');xlabel('Time');ylabel('State');legend;title('KalmanFilterSimulation');演示结果如下：5-1如何通过实验验证LMS滤波器的收敛性？解析：验证LMS算法收敛性的常用方法包括：‌绘制权值更新曲线‌：观察权值随迭代次数的变化趋势，评估收敛速度和稳定性。‌计算均方误差（MSE）‌：在收敛后，MSE应趋近于最小值，表示误差被有效抑制。‌比较不同步长μ的效果‌：通过实验验证μ对收敛速度和稳态误差的影响。‌对比LMS与RLS性能‌：在相同场景下，评估LMS的收敛速度和稳态误差。‌使用遗忘因子‌：在非平稳信号中，结合遗忘因子（如遗忘因子RLS）验证跟踪能力。5-2RLS滤波器在非平稳信号处理中的优势是什么？解析：RLS算法在非平稳信号处理中的优势包括：快速跟踪能力‌：通过遗忘因子λ<1自动调整输入权重，适应快速变化环境（如雷达目标跟踪）。高精度估计‌：稳态误差极小，尤其在平稳信号中表现优异。动态适应‌：无需预知输入统计特性，实时更新权值向量。5-3LMS与RLS滤波器的主要区别是什么？解析：LMS算法与RLS算法的主要区别包括：‌收敛速度‌：RLS收敛速度快（通常在几个采样周期内），LMS收敛慢。‌稳态误差‌：RLS稳态误差极小，LMS稳态误差较大。‌计算复杂度‌：RLS计算复杂度为O(M²)，LMS计算复杂度为O(M)。‌适应非平稳信号‌：RLS通过遗忘因子λ<1自动调整输入权重，LMS需结合遗忘因子或自适应步长提升跟踪能力。‌数值稳定性‌：RLS需使用平方根RLS（如Cholesky分解）保持数值稳定性，LMS无需矩阵运算。5-4已知自适应滤波器输入信号的自相关矩阵的最大特征值为，目标信号的均方值为。（1）求最陡下降法的步长参数的取值范围。（2）若步长参数，判断算法是否稳定。解析：（1）最陡下降法的收敛条件为：代入，得：因此，的取值范围为。（2）当时，等于收敛条件的上限，算法可能不稳定。需进一步分析具体信号特性。5-5给定自适应滤波器输入信号，目标信号，初始权重，步长参数。（1）设计滤波器阶数N=2，计算前5次迭代的误差。（2）画出滤波器的学习曲线（误差随迭代次数变化）。解析：（1）前5次迭代误差（计算过程略）：（2）学习曲线图（示例）：迭代次数|误差e(n)0|1.001|0.982|0.963|0.944|0.925-6假设自适应滤波器输入信号为自然模式（即为单位矩阵），目标信号的均方值为。（1）求最陡下降法的瞬态误差梯度。（2）若步长参数，计算第1次迭代后的权重。解析：（1）瞬态误差梯度：其中为输入信号，为目标信号共轭。（2）第1次迭代（假设初始权重）：梯度：权重更新：5-7已知自适应滤波器采用LMS算法，滤波器长度N=2，期望信号d(n)、输入信号x(n)的序列如下表，步长因子μ=0.1，初始权向量w(0)=[0,0]T。输入信号x(n)与期望信号d(n)序列n012x(n)123d(n)358按照题目要求计算n=0,1,2时刻的权向量w(n)、输出信号y(n)和误差信号e(n)。解析：Step1.n=0时刻输入向量：x(0)=[x(0),x(−1)]T=[1,0]T输出信号：y(0)=wT(0)x(0)=0×1+0×0=0误差信号：e(0)=d(0)−y(0)=3−0=3权向量更新：w(1)​=w(0)+2μe(0)x(0)=[0,0]T+2×0.1×3×[1,0]T=[0.6,0]T​Step2.n=1时刻输入向量：x(1)=[x(1),x(0)]T=[2,1]T输出信号：y(1)=wT(1)x(1)=0.6×2+0×1=1.2误差信号：e(1)=d(1)−y(1)=5−1.2=3.8权向量更新：w(2)=w(1)+2μe(1)x(1)=[0.6,0]T+2×0.1×3.8×[2,1]T=[0.6+1.52,0+0.76]T=[2.12,0.76]T​Step3.n=2时刻输入向量：x(2)=[x(2),x(1)]T=[3,2]T输出信号：y(2)=wT(2)x(2)=2.12×3+0.76×2=7.98误差信号：e(2)=d(2)−y(2)=8−7.98=0.02权向量更新：w(3)=w(2)+2μe(2)x(2)=[2.12,0.76]T+2×0.1×0.02×[3,2]T=[2.12+0.012,0.76+0.008]T=[2.132,0.768]T​最终结果：nw(n)x(n)y(n)e(n)w(n+1)0[0,0]T[1,0]T03[0.6,0]T1[0.6,0]T[2,1]T1.23.8[2.12,0.76]T2[2.12,0.76]T[3,2]T7.980.02[2.132,0.768]T随着迭代次数增加，误差信号e(n)从3快速减小到0.02，表明LMS算法的权向量在不断逼近最优值，自适应滤波效果逐步提升。5-8已知自适应滤波器采用RLS算法，滤波器长度N=2，遗忘因子λ=1（即最小二乘算法，无遗忘），初始权向量w(0)=[0,0]T，初始逆相关矩阵P(0)=δ-1输入信号x(n)与期望信号d(n)序列如下：n012x(n)123d(n)358按要求计算n=0,1,2时刻的权向量w(n)、输出信号y(n)、误差信号e(n)。解析：Step1.n=0时刻输入向量：x(0)=[x(0),x(−1)]T=[1,0]T输出信号：y(0)=wT(−1)x(0)，因w(−1)=w(0)=[0,0]T，故y(0)=0×1+0×0=0先验误差：e(0)=d(0)−y(0)=3−0=3增益向量：k(0)=逆相关矩阵更新：P权向量更新：w(0)=w(Step2.n=1时刻输入向量：x(1)=[x(1),x(0)]T=[2,1]T输出信号：y(1)=wT(0)x(1)=1.5×2+0×1=3先验误差：e(1)=d(1)−y(1)=5−3=2增益向量：xk逆相关矩阵更新：P权向量更新：wStep3.n=2时刻输入向量：x(2)=[x(2),x(1)]T=[3,2]T输出信号：y(2)=wT(1)x(2)=2×3+0.5×2=7先验误差：e(2)=d(2)−y(2)=8−7=1增益向量：xT2权向量更新：w最终结果：nw(n−1)x(n)y(n)e(n)w(n)0[0,0]T[1,0]T03[1.5,0]T1[1.5,0]T[2,1]T32[2,0.5]T2[2,0.5]T[3,2]T71[2.125,0.875]TRLS算法通过递归更新权向量，收敛速度远快于LMS算法，从结果可以看到，3次迭代后误差已降至1，权向量快速逼近最优解；且遗忘因子λ=1时，RLS等价于批处理最小二乘算法。5-9用MATLAB中的dsp.LMSFilter和dsp.RLSFilter函数分别实现LMS和RLS自适应滤波器的仿真演示，并分析误差结果。解析：（1）LMS的MATLAB代码示例如下：%生成测试信号t=0:0.001:1;d=sin(2*pi*50*t)+0.5*randn(size(t));%含噪声的期望信号x=sin(2*pi*50*t);%参考输入信号%创建LMS滤波器lmsFilter=dsp.LMSFilter('Length',32,'StepSize',0.01);[y,e,w]=lmsFilter(x',d');%绘制结果figure;subplot(2,1,1);plot(d);title('原始信号+噪声');subplot(2,1,2);plot(e);title('滤波后信号');仿真结果演示如下：（2）RLS的MATLAB代码示例如下：%生成测试信号（同上）d=sin(2*pi*50*t)+0.5*randn(size(t));x=sin(2*pi*50*t);%创建RLS滤波器rlsFilter=dsp.RLSFilter('Length',32,'ForgettingFactor',0.98);[y,e]=rlsFilter(x',d');%绘制误差曲线figure;plot(abs(e));title('RLS算法误差收敛曲线');仿真结果演示如下：5-10试写一个简单的MATLAB代码示例，演示如何实现LMS滤波器算法，绘出输入/输出信号对比及均方误差曲线，并分析仿真结果。解析：MATLAB代码示例如下：%%LMS算法MATLAB实现clear;clc;closeall;%==========参数设置==========g=100;%统计仿真次数N=1024;%输入信号抽样点数k=128;%滤波器阶数u=0.0002;%步长因子%==========发射端处理==========forq=1:gt=1:N;a=1;s=a*sin(0.05*pi*t);%输入单频信号xn=awgn(s,5);%加入高斯白噪声（SNR=5dB）%==========初始化==========y=zeros(1,N);%输出信号y(1:k)=xn(1:k);%初始输出为输入前k点w=zeros(1,k);%抽头系数初值e=zeros(1,N);%误差信号%==========LMS迭代滤波==========fori=(k+1):NXN=xn((i-k+1):i);%输入向量y(i)=w*XN';%滤波器输出e(i)=s(i)-y(i);%误差信号w=w+u*e(i)*XN;%抽头系数更新end%==========存储误差==========pp(q,:)=(e(k+1:N)).^2;%存储每次仿真的误差平方end%==========统计平均误差==========forb=1:(N-k)bi(b)=sum(pp(:,b))/g;%误差统计平均end%==========绘图==========figure;subplot(3,1,1);plot(t,real(s));title('信号s时域波形');xlabel('n');ylabel('s');axis([0N-a-1a+1]);subplot(3,1,2);plot(t,real(xn));title('信号s加噪声后的时域波形');xlabel('n');ylabel('s+噪声');subplot(3,1,3);plot(t,real(y));title('自适应滤波后的输出时域波形');xlabel('n');ylabel('y');figure;t=1:(N-k);plot(t,bi,'r');title('算法收敛曲线');xlabel('迭代次数');ylabel('均方误差');holdon;结果演示：输入/输出对比：LMS算法可有效去除信号中的噪声，输出信号收敛到输入信号。收敛曲线：均方误差随迭代次数增加逐渐降低，最终趋于稳态。步长影响：步长μ越大，收敛速度越快，但稳态误差越高；反之则收敛慢但误差低。6-1已知实信号x(n)的频谱如图6.10.1所示，试画出该信号2倍抽取后的频谱X图6.10.1题6-1图解：xX也可以写为：X6-2已知实信号x(n)图6.10.2题6-2图解：用公式分段表示（以ω>Xω<X6-3已知实信号x(n)图6.10.3题6-3图解：设X验证:负频：X验证：ω=-0.5π得1，6-4试求图6.10.4所示多速率系统的输入x(n)图6.10.4题6-4图解析：各模块均为整数倍插值或抽取，无滤波器。第一步:v1即：v1第二步:vv2代入v1v冲激函数非零当且仅当：3n由于k必须是整数，所以n必须是7的倍数。令n=7lv即v2n仅在n=第v3代入v2的特性：v2t仅在t=7时满足：1.2.n7=m得n=v其它n处v3n=0v第y已知v3t仅在t=1.n=2.n3=得n=y其他n处y(n)=0。因此：y等价地，可写作：y因为当n=147l时，n6-5已知实信号x(n)的频谱如图6.10.5所示，试画出该信号3倍抽取后的频谱。图6.10.5题6-5图解：XX即：X是一个V形，在ω=0时为0,在6-6抽取器如图6.10.6(a)所示，已知某模拟信号的频谱如图6.10.6(b)所示,其最高频率为，若对该信号以。进行采样,形成，然后通过图6.10.6(a)所示的抽取器,得到抽取后的信号,试分析并分别画出以下波形：原序列的频谱(以为变量)，原序列的频谱(以为变量)；抽取序列的频谱(以为变量)；抽取序列的频谱(以=为变量)。图6.10.6题6-6图（1）由于XajΩ是三角，支撑在-ΩXejΩT(以Ω为变量)=并在Ω>Ωh（2）频谱为：X并在ωh<ω≤（3）X并在(4)6-7在图6.10.7所示系统中，已知分别是理想的实系数的低通、带通和高通滤波器，通带为1，阻带为零，通带频率分别是，已知的频率响应如图6.10.7(b)所示，分别画出的幅频响应。图6.10.7题6-7图解：X负频率ω<因为Xejω是以2π为周期的，可以画出主值区间·H·H·H1H0通带ω≤πY(2)y1H1通带π/3≤ω≤2π/3Y3H2通带2π/3≤ω≤π(正频)及-Y6-8已知序列是由图6.10.8(a)所示的系统得到，序列是由图6.10.8(b)所示的系统得到，图题6.10.8(c)是模拟滤波器的频率特性。现希望用数字域方法直接从得到，试给出具体实现方法的框图。实现中用到数字滤波器时，给出具体指标要求。图6.10.8题6-8图解：·输入：xa·模拟低通滤波器HH其中f=Ω·采样间隔：T1·采样率：fs·输出序列：x·输入：x·模拟低通滤波器HH·采样间隔：T2·采样率：fs·输出序列：x已知xn希望得到x2采样率变换比：f这是分数倍采样率下降,可采用整数倍内插-抽取法实现。（1）内插-抽取方案采用L=xn(10kHz)→↑3(2)低通滤波器Hej·工作在内插后采样率fs·模拟带宽限制：fmax=3000·对应数字截止频率:ω·通带:0≤ω≤·阻带：至少从0.2π到π足够衰减，以消除：1.内插产生的镜像频谱2.防止后续抽取时混叠·滤波器类型：可用FIR或IIR，线性相位推荐FIR。6-9数字录音带（DAT）驱动器的采样率为48kHz，而光盘(CD)播放机则以44.1kHz的采样率工作。为了直接把声音从CD录制到DAT，需要把采样率从44.1kHz转换到48kHz。为此，考虑采用图6.10.9所示的分数倍采样率变换系统。图6.10.9题6-9图求L和M的最小可能值，以及低通滤波器H(ejw)的解：采样率转换比：f系统采样率变化关系：M为使L与M互素的最小正整数解，取L为实现分数倍变换160147x其中Hejωf滤波器的两个作用：抗镜像滤波和抗混叠滤波采样率：7.056MHz通带截止频率fp保守设计取：CD信号带宽为22.05kHz(Nyquist频率),通常音频上限取20kHz。f阻带起始频率f输出Nyquist频率为24kHz,为避免抽取混叠，阻带应在24kHz之前开始。可设：f数字频率对应：ωω过渡带非常窄，要求滤波器具有很高的阶数和陡峭的滚降。衰减要求：通带纹波≤0.1dB6-10某插值系统如图6.10.10(a)所示，系统中是图6.10.10(a)右边所示的5点序列。图6.10.10(b)是图6.10.10(a)的等效实现，其中的长度不超过3点()。试求：对任意给定的，正确选择与的关系，使两系统等效，即输出=。图6.10.10题6-10图解：

根据h(n)长度5确定多相分量设hn={h按L=3p即p0p即p1p2m=即p2=为了结构对称，可能h1n,h2系统(b)输入xn·上分支：xn→h1·中分支：xn→h2·下分支：xn→h3然后w2n延迟1,w3系统(a)的输出：y而y=要求对任意xn有y令rh这里r=对hnnnnnnn所以：hl(n)h2h37-1现有一低通滤波器，其频率响应为：H(ω)且满足G(ω)=-exp(-jω)H*解：因为H(ω)是实值的(0或1),所以HG(ω)=-H(ω)的支撑集（非零区间）为：|ω|<那么H(ω+π)的非零区间为：|ω+π|<π即：--所以：H(ω+π)=于是：G(ω)=直接用逆傅里叶变换求小波ψ(t)ψ(t)=代入G(ω)=-e-jω⋅1ψ(t)=ψ(t)=-ψ(t)=-ψ(t)=-7-2现利用Harr小波对一信号观测序列{3,5,9,11,13,15,17,19}进行三级多分辨率分析，请给出每一级的变换结果。解：低频近似：a高频细节：d第1级分解：3ad(9,11):ad(13,15):a2d(17,19):a3d所以：Level1近似（低频）aLevel1细节（高频）d第2级分解：现在对a1(4,10):a0d(14,18):a1d所以：Level2近似aLevel2细节d第3级分解：现在对a2(7，16):a0d所以：Level3近似aLevel3细节d7-3设信号为x(t)=Acos(ω0t)解：WTMorlet小波（复值）通常取为：ψ(t)=C则ψ则WT又因为cos则WT忽略常数中π1/4WT7-4设φ(t),ψ(t)是Harr尺度和小波函数，信号x(t)存在于V0子空间，并定义如下：x(t)将x(t)分解到子空间W1、W2、V2,分别求各小波和尺度系数。解：Haar尺度函数：φ(t)=Haar小波函数：ψVW信号x(t)∈V0，所以可以用c其中V因此V我们要将使用非归一化的平均/差形式ad其中第一层分解（V0a所以第二层分解（V1对a1ad所以可以整体写出：x(t)=0.5⋅其中7-5设信号x(t)=e-t2/10(sin2t+2cos4t+0.5sin(t)sin(50t)),用间隔T=1/28，从t用MATLAB的DWT函数，将an0分解到W1，…，W6，V用IDWT函数，进行合成实验，重构an0；假设在an解：分解并画图步骤：1.生成信号x(t)的采样。2.采样间隔T=1/28,从3.用dwt或wavedec进行多级分解到第6层，得到W14.分别画出各层小波系数和最后一层尺度系数。MATLAB代码：%参数设置N=256;T=1/256;%1/2^8t=(0:N-1)*T;%生成信号x=exp(-t.^2/10).*(sin(2*t)+2*cos(4*t)+0.5*sin(t).*sin(50*t));%设定小波类型，比如db4wavelet='db4';level=6;%多级分解[C,L]=wavedec(x,level,wavelet);%提取各层系数approx_coef=appcoef(C,L,wavelet,level);%V6系数detail_coefs=cell(1,level);fork=1:leveldetail_coefs{k}=detcoef(C,L,k);%W_k系数end%画图：近似系数（V6）figure;subplot(level+1,1,1);plot(approx_coef);title(['ApproximationCoefficientsV'num2str(level)]);xlabel('Position');ylabel('Amplitude');gridon;%画图：细节系数W1到W6fork=1:levelsubplot(level+1,1,k+1);plot(detail_coefs{k});title(['DetailCoefficientsW'num2str(k)]);xlabel('Position');ylabel('Amplitude');gridon;endsgtitle('WaveletDecompositionCoefficients');分解画图：按上述代码画出V6第（2）问：合成重构步骤：直接使用

waverec（或

idwt

逐层重构）对上面分解的系数进行重构，验证是否与原信号一致。MATLAB代码：%重构信号x_recon=waverec(C,L,wavelet);%计算重构误差err=max(abs(x-x_recon));disp(['最大重构误差:',num2str(err)]);%画图比较figure;subplot(2,1,1);plot(t,x,'b','LineWidth',1.5);title('OriginalSignal');xlabel('t');ylabel('x(t)');gridon;subplot(2,1,2);plot(t,x_recon,'r--','LineWidth',1.5);title('ReconstructedSignal');xlabel('t');ylabel('x_{recon}(t)');gridon;sgtitle('OriginalvsReconstructedSignal');重构实验:用waverec重构,误差极小(<10第（3）问：小波去噪去噪算法设计步骤：1.对含噪信号进行小波分解到第6层。2.对每一层细节系数进行阈值处理(硬阈値或软阈值)。3.阈值选择：通用阈值(UniversalThreshold)σ2log⁡N用第一层细节系数的中位数绝对值估计：σ4.用处理后的系数重构信号。MATLAB代码：%添加噪声rng(0);noise=randn(1,N);x_noisy=x+noise;%分解含噪信号[Cn,Ln]=wavedec(x_noisy,level,wavelet);%估计噪声标准差（从第一层细节系数）detail1=detcoef(Cn,Ln,1);sigma=median(abs(detail1))/0.6745;%通用阈值universal_thr=sigma*sqrt(2*log(N));%对各层细节系数进行软阈值处理thresholded_C=Cn;%复制系数向量start_idx=Ln(1)+1;%跳过近似系数部分的位置fork=1:levellength_detail=Ln(level+2-k);%细节系数长度detail=thresholded_C(start_idx:start_idx+length_detail-1);%软阈值detail=sign(detail).*max(abs(detail)-universal_thr,0);thresholded_C(start_idx:start_idx+length_detail-1)=detail;start_idx=start_idx+length_detail;end%重构去噪信号x_denoised=waverec(thresholded_C,Ln,wavelet);%画图比较figure;subplot(3,1,1);plot(t,x,'b');title('OriginalCleanSignal');gridon;subplot(3,1,2);plot(t,x_noisy,'m');title('NoisySignal(GaussianWhiteNoise,var=1)');gridon;subplot(3,1,3);plot(t,x_denoised,'r','LineWidth',1.5);title('DenoisedSignal(WaveletSoftThresholding)');gridon;%计算信噪比改善SNR_noisy=20*log10(norm(x)/norm(x_noisy-x));SNR_denoised=20*log10(norm(x)/norm(x_denoised-x));disp(['SNR_noisy:',num2str(SNR_noisy),'dB']);disp(['SNR_denoised:',num2str(SNR_denoised),'dB']);去噪算法：采用小波软阈值去噪，阈值用通用阈值，噪声标准差由第一层细节系数中位数估计。去噪后信号SNR提高，视觉效果改善。7-6对应于Daub4小波的低通与高通分解滤波器如下LO_D=-0.12940.22410.83650.4830HO_D=-0.48300.8365-0.2241-0.1294证明这两个滤波器自身正交。证明：设h[n]为低通滤波器，支撑集为n=0,1,2,3。设g[n]为高通滤波器，支撑集也为n=0,1,2,3。要证明的是分析滤波器h与g在偶数平移上正交，即：∀由于滤波器有限长，求和只需考虑n使h[n]和g[n-2k]都非零的项。先验证是否满足Daubechies正交小波的镜像对称公式:g[n]=(-1逐一验证：∙n=0:(-1∙n=1:(-1∙n=2:(-1∙n=3:(-1成立。因此g是h的交替翻转形式（带一个全局符号-1对于所有n实际等价于(-1)^{n+1})h[n]只在n=0,1,2,3非零。g所以重叠的条件是：0≤n≤3,2k≤n≤2k+3.这意味着：max(0,2k)≤n≤min(3,2k+3).K可以取的使重叠非空的值有限：·k=0:区间[0,3]与[0,3]重叠长度为4·k=1:区间[0,3]与[2,5]重叠：n=2,3·k=-1:区间[0,3]与[-2,1]重叠：n=0,1·其他k无重叠，和为0显然。因此只需要验证k=-1,0,1&&(&(对于所有整数k,都有n因此这两个滤波器h与g满足正交条件，即它们自身正交(分析滤波器之间的偶数平移正交性)。7-7利用MATLAB提供的chirp函数产生一个chirp信号，其频率从50Hz变化到1000Hz，持续时间为2s，信号的采样频率为2500Hz。利用Morlet小波对信号进行连续小波变换，显示其时间-频率分布图；采样频率不变，当chirp信号的频率在2s内从50Hz变化到2000Hz，再利用Morlet小波对信号进行连续小波变换，诠释出现混叠的原因。解：（1）MATLAB代码%%参数设置fs=2500;%采样频率2500HzT=2;%信号持续时间2st=0:1/fs:T-1/fs;%时间向量N=length(t);%采样点数=5000%%1.生成第一个chirp信号(50Hz→1000Hz)chirp1=chirp(t,50,T,1000,'linear');%%2.Morlet小波连续小波变换figure;[cwt_coeffs1,freq1]=cwt(chirp1,'amor',fs);%'amor'表示复Morlet小波%绘制时频分布图imagesc(t,freq1,abs(cwt_coeffs1));axisxy;colormap(jet);colorbar;xlabel('时间(s)');ylabel('频率(Hz)');title('Chirp信号时频图(50Hz→1000Hz)');set(gca,'YScale','log');%对数频率轴便于观察结果说明：第一个信号的时频图显示频率从50Hz线性增加到1000Hz，频率轨迹清晰连续，没有异常现象。因为信号最高频率（1000Hz）小于Nyquist频率（1250Hz），满足采样定理，无混叠。（2）MATLAB代码%%3.生成第二个chirp信号(50Hz→2000Hz)chirp2=chirp(t,50,T,2000,'linear');%%4.Morlet小波连续小波变换figure;[cwt_coeffs2,freq2]=cwt(chirp2,'amor',fs);%绘制时频分布图imagesc(t,freq2,abs(cwt_coeffs2));axisxy;colormap(jet);colorbar;xlabel('时间(s)');ylabel('频率(Hz)');title('Chirp信号时频图(50Hz→2000Hz)-出现混叠');set(gca,'YScale','log');结果说明：第二个信号的时频图出现异常：频率轨迹在约1.5秒后发生畸变，高频部分（>1250Hz）的能量出现在低频区域，形成虚假的频率成分。本实验中：·采样频率f·Nyquist频率f·第二个信号的最高频率f因此违反了采样定理，导致混叠。混叠是采样过程中违反Nyquist定理的直接后果，时频分析方法（包括小波变换）只能显示采样后信号的特征，无法恢复混叠前的高频信息。7-8利用Daub4小波演示一个单位脉冲信号通过一级离散小波变换的全过程，并利用MATLAB提供的dwt和idwt函数验证其结果；分析离散小波变换后的数据，以及重建后的数据。（可利用MATLAB函数wfilters获得Daub4小波对应的滤波器系数，即[LoDHiDLoRHiR]=wfilters('db2')）解：MATLAB代码%%1.获取Daub4(db2)滤波器系数[LoD,HiD,LoR,HiR]=wfilters('db2');disp('Daub4滤波器系数:');disp('低通分解LoD:');disp(LoD');disp('高通分解HiD:');disp(HiD');disp('低通重构LoR:');disp(LoR');disp('高通重构HiR:');disp(HiR');%%2.生成单位脉冲信号x=zeros(1,10);x(5)=1;%δ[n-4](MATLAB索引从1开始)disp('原始信号x:');disp(x);%%3.使用dwt函数进行一级分解[cA1,cD1]=dwt(x,LoD,HiD);disp('近似系数cA1(dwt输出):');disp(cA1);disp('细节系数cD1(dwt输出):');disp(cD1);%%4.手工计算验证%卷积后下采样conv_low=conv(x,LoD,'full');conv_high=conv(x,HiD,'full');cA1_manual=conv_low(2:2:end);%MATLABdwt从索引2开始取（保持相位）cD1_manual=conv_high(2:2:end);disp('手工计算cA1:');disp(cA1_manual);disp('手工计算cD1:');disp(cD1_manual);%检查与dwt是否一致fprintf('cA1匹配误差:%e\n',norm(cA1-cA1_manual));fprintf('cD1匹配误差:%e\n',norm(cD1-cD1_manual));%%5.使用idwt重构x_rec=idwt(cA1,cD1,LoR,HiR);disp('重建信号x_rec:');disp(x_rec);disp('重建误差:');disp(norm(x-x_rec));运行结果：Daub4滤波器系数:低通分解LoD:[-0.1294;0.2241;0.8365;0.4830]高通分解HiD:[-0.4830;0.8365;-0.2241;-0.1294]低通重构LoR:[0.4830;0.8365;-0.2241;0.1294]高通重构HiR:[-0.1294;-0.2241;0.8365;-0.4830]原始信号x:[0000100000]近似系数cA1(dwt输出):[0.22410.4830000]细节系数cD1(dwt输出):[0.8365-0.1294000]手工计算cA1:[0.22410.483000000]手工计算cD1:[0.8365-0.129400000]cA1匹配误差:0.000000e+00cD1匹配误差:0.000000e+00重建信号x_rec:[0000100000]重建误差:1.5701e-167-9分别采用Daub4小波和Haar小波，利用MATLAB对下列定义的信号进行离散小波变换，显示其一级近似系数和细节系数。fs=10000;T=1/fst=0:T:4095*T;fn=cos(2*pi*6*t);fn(1250:1255)=fn(1249);该信号是一个频率为6Hz的余弦信号，其在样点1250~1255之间有个毛刺干扰。解：MATLAB代码%%1.信号生成fs=10000;%采样频率10000HzT=1/fs;%采样间隔t=0:T:4095*T;%时间向量，4096个点fn=cos(2*pi*6*t);%6Hz余弦信号fn(1250:1255)=fn(1249);%毛刺干扰（样点1250~1255）%%2.获取小波滤波器系数%Daub4小波（db2）[LoD_d4,HiD_d4,~,~]=wfilters('db2');%Haar小波（db1）[LoD_haar,HiD_haar,~,~]=wfilters('db1');%%3.一级离散小波变换%使用Daub4小波[cA1_d4,cD1_d4]=dwt(fn,LoD_d4,HiD_d4);%使用Haar小波[cA1_haar,cD1_haar]=dwt(fn,LoD_haar,HiD_haar);%%4.显示结果figure('Position',[100,100,1200,800]);%(1)原始信号subplot(4,2,1);plot(t,fn);xlabel('时间(s)');ylabel('幅度');title('原始信号(6Hz余弦+毛刺)');gridon;xlim([0,0.5]);%显示前0.5秒%标记毛刺位置holdon;plot(t(1249:1256),fn(1249:1256),'r','LineWidth',2);holdoff;%(2)Daub4近似系数subplot(4,2,3);plot(1:length(cA1_d4),cA1_d4);xlabel('样点索引');ylabel('幅度');title('Daub4小波-近似系数(cA1)');gridon;%(3)Daub4细节系数subplot(4,2,4);plot(1:length(cD1_d4),cD1_d4);xlabel('样点索引');ylabel('幅度');title('Daub4小波-细节系数(cD1)');gridon;%标记毛刺对应位置（下采样后）holdon;dwt_loc=ceil(1250/2);%DWT后毛刺大致位置plot(dwt_loc,cD1_d4(dwt_loc),'ro','MarkerSize',8,'LineWidth',2);holdoff;%(4)Haar近似系数subplot(4,2,5);plot(1:length(cA1_haar),cA1_haar);xlabel('样点索引');ylabel('幅度');title('Haar小波-近似系数(cA1)');gridon;%(5)Haar细节系数subplot(4,2,6);plot(1:length(cD1_haar),cD1_haar);xlabel('样点索引');ylabel('幅度');title('Haar小波-细节系数(cD1)');gridon;%标记毛刺对应位置holdon;plot(dwt_loc,cD1_haar(dwt_loc),'ro','MarkerSize',8,'LineWidth',2);holdoff;%(6)信号频谱subplot(4,2,2);N=length(fn);f=fs*(0:N/2-1)/N;spectrum=abs(fft(fn));plot(f,spectrum(1:N/2));xlabel('频率(Hz)');ylabel('幅度');title('信号频谱');gridon;xlim([0,100]);%(7)细节系数对比subplot(4,2,7:8);plot(1:length(cD1_d4),cD1_d4,'b','LineWidth',1.5);holdon;plot(1:length(cD1_haar),cD1_haar,'r','LineWidth',1.5);xlabel('样点索引');ylabel('幅度');title('细节系数对比(蓝色:Daub4,红色:Haar)');gridon;legend('Daub4','Haar');holdoff;运行结果与分析1.系数长度原始信号长度：4096点；一级DWT后系数长度：2048点（下采样2倍）2.Daub4小波结果近似系数cA1：平滑地表示6Hz余弦信号的低频成分细节系数cD1：在毛刺位置（样点625附近）有明显的脉冲响应Daub4小波长度=4，毛刺影响范围较宽毛刺能量分散到多个系数中3.Haar小波结果近似系数cA1：阶跃式表示信号的低频成分细节系数cD1：在毛刺位置有明显的脉冲响应Haar小波长度=2，毛刺响应更局部毛刺能量集中在少数几个系数中7-10基于MATLAB，利用Daub4小波对图像进行多尺度统计分析，要求：(1)利用网络资源搜集10幅质量较高的图像，要求图像分辨率≥256×256像素,保存为.jpg或.png格式;(2)对每幅图像进行4层二维小波分解，生成13个子带;(3)计算每个子带的小波系数均值和方差。解析：1）二维小波分解：对图像行和列分别进行一维DWT；每层分解产生4个子带：近似(A)、水平(H)、垂直(V)、对角(D)；4层分解共产生13个子带2）Daub4小波特性：MATLAB函数：wfilters('db2')；滤波器长度：4；正交性：满足完美重建条件统计量计算均值：反映子带系数的平均强度μ=方差：反映子带系数的离散程度σ（说明：在当前目录创建images文件夹放入10幅相同图像(.png格式））MATLAB代码%%基于Daub4小波的图像统计分析clear;clc;closeall;%%1.设置image_folder='images';wavelet_name='db2';level=4;%%2.检查文件夹if~exist(image_folder,'dir')mkdir(image_folder);fprintf('已创建文件夹:%s\n请放入10幅图像后重新运行程序。\n',image_folder);return;end%%3.获取图像文件files=dir(fullfile(image_folder,'*.jpg'));files=[files;dir(fullfile(image_folder,'*.png'))];iflength(files)<1fprintf('请在images文件夹中放入jpg或png格式的图像。\n');return;endnum_images=min(length(files),10);files=files(1:num_images);%%4.子带名称bands={'A4','H4','V4','D4','H3','V3','D3','H2','V2','D2','H1','V1','D1'};num_bands=13;%%5.初始化mean_val=zeros(num_images,num_bands);var_val=zeros(num_images,num_bands);%%6.处理图像fori=1:num_imagesfprintf('处理%d/%d:%s\n',i,num_images,files(i).name);%读取图像img=imread(fullfile(image_folder,files(i).name));ifsize(img,3)==3img=rgb2gray(img);endimg=im2double(img);%调整大小ifsize(img,1)<256||size(img,2)<256img=imresize(img,[256256]);end%小波分解[C,S]=wavedec2(img,level,wavelet_name);%提取系数idx=1;%A4A4=appcoef2(C,S,wavelet_name,level);mean_val(i,idx)=mean(A4(:));var_val(i,idx)=var(A4(:));idx=idx+1;%各层细节forlev=level:-1:1[H,V,D]=detcoef2('all',C,S,lev);mean_val(i,idx)=mean(H(:));var_val(i,idx)=var(H(:));idx=idx+1;mean_val(i,idx)=mean(V(:));var_val(i,idx)=var(V(:));idx=idx+1;mean_val(i,idx)=mean(D(:));var_val(i,idx)=var(D(:));idx=idx+1;endend%%7.结果显示fprintf('\n=