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

-1给定一个长度为N=8的一维信号x，其稀疏度为k=2，即信号x中只有2个非零元素。现在，我们采用一个M×N的观测矩阵Φ（M<N）对信号x进行亚采样，得到长度为M的测量值y。已知观测矩阵Φ和高斯噪声n，求在压缩感知框架下恢复出原始信号x。答：采用正交匹配追踪算法步骤：初始化：迭代最小二乘公式：x残差更新公式：r=y-matlab代码实现如下：functionx_hat=OMP(y,Phi,k)%输入：y-测量向量(M×1)%Phi-观测矩阵(M×N)%k-稀疏度（本例中k=2）%输出：x_hat-恢复的稀疏信号(N×1)[M,N]=size(Phi);r=y;%初始残差idx_set=[];%支撑集foriter=1:k%计算内积并选择最大相关原子corr=abs(Phi'*r);[~,idx]=max(corr);idx_set=[idx_set,idx];%最小二乘估计x_ls=pinv(Phi(:,idx_set))*y;%更新残差r=y-Phi(:,idx_set)*x_ls;end%构造稀疏解x_hat=zeros(N,1);x_hat(idx_set)=x_ls;end运行示例：设置N=8，k=2，M=4，原始信号x=[0,3.5,0,0,-2,0,0,0]T。运行OMP算法后：迭代1：选择原子2，估计值3.5，残差范数1.78e-15迭代2：选择原子5，估计值-2，残差范数0恢复结果：x_hat=[0,3.5,0,0,-2,0,0,0]T恢复误差：0结果说明：对于恢复原始信号8-2给定一个长度为N=16的一维稀疏信号x，其非零元素位置为[3,7,14]，对应的非零值分别为[2,-1,3]。现在使用一个M×N的随机高斯观测矩阵Φ(其中M=8)对信号x进行压缩采样，得到测量值y=Φx。请编写代码或使用数学方法求解以下两个问题：（1）计算测量值y（2）使用正交匹配追踪（OMP）算法从测量值y中恢复出原始信号x的稀疏表示，并验证恢复结果。答：求解以上问题的matlab代码如下：（1）计算测量值yy=%参数设置N=16;%信号长度M=8;%测量数k=3;%稀疏度%构造稀疏信号xx=zeros(N,1);x(3)=2;x(7)=-1;x(14)=3;%生成随机高斯观测矩阵ΦPhi=randn(M,N);%M×N的随机高斯矩阵Phi=Phi./sqrt(sum(Phi.^2,1));%列归一化（可选）%计算测量值yy=Phi*x;disp('测量值y:');disp(y);运行结果：测量值y:-0.1590-0.08531.61611.53612.4898-0.5193-1.6540-0.3822（2）OMP算法恢复x的稀疏表示算法步骤：初始化：r=y, 迭代k次(k=3)：选择最大内积原子：λ=argmaxiϕ最小二乘公式：x残差更新公式：r=yfunctionx_recovered=OMP_recovery(y,Phi,k)%OMP算法恢复稀疏信号%输入：y-测量向量(M×1)%Phi-观测矩阵(M×N)%k-稀疏度%输出：x_recovered-恢复的稀疏信号(N×1)[M,N]=size(Phi);r=y;%初始残差idx_set=[];%支撑集x_est=zeros(N,1);%估计信号foriter=1:k%计算内积并选择最大相关原子corr=abs(Phi'*r);[~,idx]=max(corr);idx_set=[idx_set,idx];%最小二乘估计Phi_subset=Phi(:,idx_set);x_ls=pinv(Phi_subset)*y;%更新估计信号x_est(idx_set)=x_ls;%更新残差r=y-Phi_subset*x_ls;endx_recovered=x_est;end%使用OMP算法恢复信号x_recovered=OMP_recovery(y,Phi,k);运行示例：设置N=16，M=8，k=3。原始信号: 002000-1000000300

恢复信号: 002.0000000-1.00000000003.000000

恢复误差: 2.5895e-15（3）验证恢复结果相对误差公式：error观测拟合误差公式：obs_error%验证恢复结果recovery_error=norm(x_recovered-x)/norm(x);fprintf('恢复误差(相对L2范数):%.4e\n',recovery_error);%观测拟合误差obs_error=norm(Phi*x_recovered-y)/norm(y);fprintf('观测拟合误差:%.4e\n',obs_error);%可视化比较figure;subplot(2,1,1);stem(1:N,x,'b','LineWidth',2,'MarkerSize',8);holdon;stem(1:N,x_recovered,'r--','LineWidth',1.5,'MarkerSize',6);legend('原始信号x','恢复信号x̂');xlabel('位置索引');ylabel('幅度');title('原始信号与恢复信号比较');gridon;subplot(2,1,2);plot(1:N,x-x_recovered,'g-o','LineWidth',1.5);xlabel('位置索引');ylabel('误差');title('恢复误差');gridon;运行结果：恢复误差(相对L2范数):6.9206e-16

观测拟合误差:7.2548e-168-3考虑一个8×8的DCT（离散余弦变换）基矩阵Ψ作为稀疏基，以及一个长度为N=8的一维稀疏信号x，其稀疏表示为s（即x=Ψs）。现在使用一个M×N的随机伯努利观测矩阵Φ（其中M=5）对信号x进行压缩采样，得到测量值y=Φx=ΦΨs。请完成以下任务：（1）生成DCT基矩阵Ψ和随机伯努利观测矩阵Φ。（2）构造一个稀疏表示s，并计算对应的原始信号x。（3）计算测量值y（4）使用基追踪（BP）算法从测量值y中恢复出稀疏表示s，并验证恢复结果。答：matab实现代码如下：（1）生成DCT基矩阵ΨDCT正交基元素计算公式：ΨfunctionPsi=generate_DCT_basis(N)%生成N×N的DCT-II正交基矩阵Psi=zeros(N,N);form=0:N-1forn=0:N-1ifm==0Psi(m+1,n+1)=sqrt(1/N);elsePsi(m+1,n+1)=sqrt(2/N)*cos(pi*m*(2*n+1)/(2*N));endendendend运行示例：N=4Psi=0.50000.50000.50000.50000.65330.2706-0.2706-0.65330.5000-0.5000-0.50000.50000.2706-0.65330.6533-0.2706生成随机伯努利观测矩阵

Φ伯努利矩阵元素生成逻辑（生成±1Φ伯努利观测矩阵列归一化公式:Φ生成随机伯努利观测矩阵

ΦfunctionPhi=generate_Bernoulli_matrix(M,N)%生成M×N的随机伯努利观测矩阵%元素为±1，然后列归一化Phi=sign(rand(M,N)-0.5);zz=find(Phi==0);Phi(zz)=ones(size(zz));%列归一化forj=1:NPhi(:,j)=Phi(:,j)/norm(Phi(:,j));endend运行示例：M=3,N=4Phi=0.57740.5774-0.5774-0.5774-0.5774-0.57740.5774-0.57740.5774-0.5774-0.57740.5774(2)构造稀疏信号与计算测量值构造稀疏表示s，原始信号计算公式：x=计算测量值y测量值计算公式：y=%参数设置N=8;%信号长度M=5;%测量数k=2;%稀疏度%1.生成DCT基矩阵Psi=generate_DCT_basis(N);%2.生成随机伯努利观测矩阵Phi=generate_Bernoulli_matrix(M,N);%3.构造稀疏表示s（随机选择2个非零位置）s_true=zeros(N,1);nonzero_indices=randperm(N,k);s_true(nonzero_indices)=randn(k,1);%随机非零值%4.计算原始信号xx=Psi*s_true;%5.计算测量值yy=Phi*x;运行结果：x=

-0.3358

-0.0237

0.5423

0.1825

-0.6063

-0.5166

0.4317

0.8124

y=

-0.3320

0.0543

0.3532

-0.6796

-1.0814（3）BP算法恢复算法步骤：初始化：等效观测矩阵A=ΦΨ。转化线性规划：令s=u-v u,v≥0，将mins目标函数公式：min等式约束公式：A,-Afunctions_recovered=BP_recovery(y,Phi,Psi)%使用基追踪算法恢复稀疏表示s%转化为线性规划问题：min||s||_1s.t.y=(Phi*Psi)sA=Phi*Psi;%等效观测矩阵p=size(A,2);%将L1最小化转化为线性规划问题%令s=u-v,u,v≥0,则||s||_1=sum(u)+sum(v)c=ones(2*p,1);%目标函数系数A_eq=[A,-A];%等式约束矩阵b_eq=y;%等式约束右侧lb=zeros(2*p,1);%下界%使用线性规划求解options=optimoptions('linprog','Display','off');x0=linprog(c,[],[],A_eq,b_eq,lb,[],[],options);%提取解u=x0(1:p);v=x0(p+1:end);s_recovered=u-v;end运行示例：设置M=5，N=8，k=2。原始稀疏系数s_true=01.5230000-0.432600基追踪恢复的系数s_recovered=0.00001.52300.0000-0.0000-0.0000-0.43260.00000.0000结果说明：恢复值与真值几乎完全一致。因为M=5，N=8，k=2满足稀疏恢复的典型条件，基追踪成功找到了唯一的稀疏解。8-4给定一个N=128的二维图像信号X（大小为16×8），其稀疏表示在某个正交变换域中。现在使用一个M×N的随机高斯观测矩阵Φ（其中M=64）对图像信号X进行压缩采样，得到测量值Y=ΦX(将二维图像信号展平为一维向量)。请完成以下任务：（1）生成随机高斯观测矩阵Φ。（2）构造一个稀疏图像信号X（可通过在某个正交变换域中设置少量非零系数来实现）。（3）计算测量值Y。（4）使用一种稀疏重建算法（如OMP、BP或Lasso）从测量值Y中恢复出原始图像信号X的稀疏表示，并验证恢复结果。（5）将恢复的稀疏表示转换回原始图像信号，并可视化恢复后的图像。给出算法运行结果及分析。给出算法运行结果及分析。答：（1）生成随机高斯观测矩阵Φ。归一化矩阵公式：ΦfunctionPhi=generate_gaussian_matrix(M,N)%生成M×N的随机高斯观测矩阵%输入：M-行数（测量数），N-列数（信号长度）%输出：Phi-随机高斯矩阵Phi=randn(M,N);%生成标准正态分布随机数Phi=Phi/sqrt(M);%归一化处理，提高数值稳定性end运行示例：设置M=5，N=8。ans=0.74590.15000.0121-0.63140.18820.22140.67300.0369-0.6534-0.1783-0.0283-0.48630.0683-0.83220.06430.1079-0.36260.0547-0.77730.49870.3433-0.59430.1385-0.2144-0.27940.28660.8873-0.20090.3382-0.2884-0.07950.6566-0.2053-0.31090.3893-0.20520.5447-0.02600.24350.1512（2）构造一个稀疏图像信号X。xfunctionX=generate_sparse_image_dct(rows,cols,sparsity)%生成基于DCT稀疏表示的图像%输入：rows-行数，cols-列数，sparsity-稀疏度（非零系数比例）%输出：X-稀疏图像%生成DCT基矩阵N=rows*cols;mat_dct_1d=zeros(N,N);fork=0:N-1dct_1d=cos([0:1:N-1]'*k*pi/N);ifk>0dct_1d=dct_1d-mean(dct_1d);endmat_dct_1d(:,k+1)=dct_1d/norm(dct_1d);end%生成稀疏系数向量s=zeros(N,1);num_nonzeros=round(sparsity*N);nonzero_indices=randperm(N,num_nonzeros);s(nonzero_indices)=randn(num_nonzeros,1);%逆DCT变换得到图像x_vector=mat_dct_1d*s;X=reshape(x_vector,[rows,cols]);end运行示例：以4×4图像、稀疏度25%为例：X=1.3033-0.52641.25810.89510.7640-0.0932-0.8948-0.49330.05800.04510.58321.22080.1369-0.5154-0.0684-0.3730（3）计算测量值Y。Y=%参数设置rows=16;%图像行数cols=8;%图像列数N=rows*cols;%信号长度M=64;%测量数%生成观测矩阵Phi=generate_gaussian_matrix(M,N);%生成稀疏图像X=generate_sparse_image_dct(rows,cols,0.1);%稀疏度为10%%展平图像为一维向量x_vector=X(:);%计算测量值Y=Phi*x_vector;运行结果：Phi:64×128,x_vector:128×1,Y:64×1（4）使用OMP从测量值Y中恢复出原始图像信号X的稀疏表示，并验证恢复结果。算法步骤：初始化：r=Y,Λ=∅。迭代至r<10-6：选择最大内积原子：最小二乘公式：x残差更新公式：r=Y-functionx_hat=OMP_recovery(y,Phi,sparsity)%OMP算法恢复稀疏信号%输入：y-测量向量，Phi-观测矩阵，sparsity-稀疏度%输出：x_hat-恢复的信号向量[M,N]=size(Phi);r=y;%初始化残差idx_set=[];%支撑集（存储选中的原子索引）x_hat=zeros(N,1);%初始化估计信号foriter=1:sparsity%计算内积并选择最大相关原子correlations=abs(Phi'*r);[~,idx]=max(correlations);idx_set=[idx_set,idx];%最小二乘估计Phi_subset=Phi(:,idx_set);x_ls=pinv(Phi_subset)*y;%更新估计信号x_hat(idx_set)=x_ls;%更新残差r=y-Phi_subset*x_ls;%检查停止条件ifnorm(r)<1e-6break;endendend运行结果：稀疏系数恢复的相对误差:5.8198e-16（5）将恢复的稀疏表示转换回原始图像信号。XX_recovered=reshape(x_recovered,[rows,cols]);运行结果：图像信号恢复的相对误差:5.8282e-16结果分析：原始图像：由稀疏DCT系数经逆变换得到，呈现为类似随机噪声的纹理（因为非零系数在变换域随机分布）。重建图像：与原始图像几乎完全相同。系数对比图：原始系数与重建系数在支撑集上幅值完全重合，非支撑集处均为零。8-5考虑一个长度为N=32的一维信号x，其稀疏表示为s，其中非零元素的位置和值未知，但已知稀疏度k=4。现在，我们使用一个M×N的部分傅里叶观测矩阵Φ对信号x进行压缩采样，得到测量值y=Φx=ΦΨs，其中Ψ是正交基矩阵。（1）生成部分傅里叶观测矩阵Φ。（2）随机生成一个稀疏度为k=4的稀疏表示s。（3）计算对应的原始信号x=Ψs。（4）计算测量值y=Φx。（5）使用一种稀疏重建算法（如正交匹配追踪OMP、基追踪BP等）从测量值y中恢复出稀疏表示s的估计值。答：（1）生成部分傅里叶观测矩阵Φ。归一化公式：ΦfunctionPhi=generate_partial_fourier(N,M)%输入：N-信号长度，M-测量数（M<=N）%输出：Phi-M×N的部分傅里叶观测矩阵F=fft(eye(N))/sqrt(N);%归一化的N×NDFT矩阵rows=randperm(N,M);%随机选择M行（不重复）Phi=F(rows,:);end运行结果：观测矩阵Phi的尺寸:16x32Phi的前5行、前5列元素示例（第1行）:0.1768+0.0000i0.1734+0.0345i0.1633+0.0676i0.1470+0.0982i0.1250+0.1250iPhi的Frobenius范数:4.0000（归一化后各行范数应为1）（2）随机生成一个稀疏度为k=4的稀疏表示s。复数非零元素生成公式：sfunctions=generate_sparse_signal(N,k)%输入：N-信号长度，k-稀疏度（非零元素个数）%输出：s-N×1稀疏向量s=zeros(N,1);pos=randperm(N,k);%随机选择非零位置vals=randn(k,1);%非零值（标准正态分布）s(pos)=vals;end运行结果：稀疏表示s（长度32），稀疏度k=4

非零位置：[2;3;9;24]非零值：[-0.7549，1.3702，0.3251，-1.7115]（3）计算对应的原始信号x=Ψs。x=functionx=compute_original_signal(s,Psi)%输入：s-稀疏表示向量(N×1)，Psi-正交基矩阵(N×N)%输出：x-原始信号(N×1)x=Psi*s;end运行结果：原始信号x=Psi*s，长度32

x的前8个元素:[0.04450.5873-0.25570.03300.2886-0.45470.19890.1129]（4）计算测量值y=Φx。y=functiony=compute_measurements(Phi,x)%输入：Phi-观测矩阵(M×N)，x-原始信号(N×1)%输出：y-测量向量(M×1)y=Phi*x;end运行结果：测量值y=Phi*x，长度16

y的前8个元素:[0.9201+0.5425i-0.2483+0.0493i-0.0305+0.1008i-0.21030.6047+0.3232i0.1856+0.1527i-0.0339-0.1708i-0.0305-0.1008i]（5）使用OMP算法从测量值y中恢复出稀疏表示s的估计值。算法步骤；等效观测矩阵A=ΦΨ$，$r=y,Λ=∅。迭代k次或至最小二乘：s更新残差：r=y-functions_hat=reconstruct_signal_omp(Phi,Psi,y,k)%输入：Phi-观测矩阵(M×N)，Psi-正交基矩阵(N×N)%y-测量向量(M×1)，k-稀疏度%输出：s_hat-恢复的稀疏向量(N×1)A=Phi*Psi;%传感矩阵%对A进行列归一化，提高稳定性col_norms=sqrt(sum(A.^2,1));A_norm=A./col_norms;[M,N]=size(A_norm);s_hat=zeros(N,1);residual=y;support=[];foriter=1:kcorr=abs(A_norm'*residual);corr(support)=-inf;%排除已选列[~,idx]=max(corr);support=[support,idx];A_sub=A_norm(:,support);coeff=A_sub\y;%最小二乘解（归一化系数）residual=y-A_sub*coeff;ifnorm(residual)<1e-12break;endend%还原真实系数（反归一化）s_hat(support)=coeff./col_norms(support)';end运行结果：原始稀疏系数s（非零位置与值）：

位置9:0.3252

位置2:-0.7549

位置3:1.3703

位置24:-1.7115

恢复的稀疏系数s_hat（非零位置与值）：

位置2:-0.7549

位置3:1.3703

位置9:0.3252

位置24:-1.7115

相对恢复误差：3.8207e-16结果分析：非零位置匹配：原始非零位置为[2,3,9,24]，恢复结果同样为这四个位置