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

7-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,14,18]再做(4,10):a0d(14,18):a1d所以：Level2近似aLevel2细节d第3级分解：现在对a2=[7，(7，16):a0d所以：Level3近似aLevel3细节d7-3设信号为x(t)=Acos(ω0t),设小波母函数为解析函数Morlet小波,求其CWT的表达式WT(解：WTMorlet小波（复值）通常取为：ψ(t)=C则ψ则WT又因为cos则WT忽略常数中π1/4等,因MorletWT7-4设φ(t),ψ(t)是Harr尺度和小波函数，信号x(t)存在于V0子空间，并定义如下：x(t)将x(t)分解到子空间W1、W2、V2,分别求各小波和尺度系数。解：Haar尺度函数：φ(t)=Haar小波函数：ψVW信号x(t)∈V0，所以可以用φ0,k(t)=φ(t-k)展开，系数就是区间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=0开始采样256用MATLAB的DWT函数，将an0分解到W1，…，W6，V6用IDWT函数，进行合成实验，重构an0；假设在an0中混入了方差为解：分解并画图步骤：1.生成信号x(t)的采样。2.采样间隔T=1/28,从t=03.用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,W6,…,W1的系数图，可以看到高频部分(W1,W2),主要对应50Hz第（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-13第（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,