【现代信号处理第六次作业】
现代信号处理第六次作业
1. 圆周卷积的定义如下
x 1 ( n ) ⊗ x 2 ( n ) = ∑ i = 0 N − 1 x 1 ( i ) x 2 ( ( n − i ) ) N x_1(n)\ x_2(n)=\sum_{i=0}^{N-1}{x_1(i)x_2((n-i)})_N x1(n)⊗x2(n)=i=0∑N−1x1(i)x2((n−i))N
循环卷积所得序列的长度还是N，长度不变。
其中 表示对 进行圆周右移m运算，如长度为6的序列 ，其图形及圆周移位的图形如下图所示：
（a） （b） （c）
长度为N的序列 通过长度为L（ L y ( n ) = x 1 ( n ) ⊗ h 1 ( n ) = ∑ i = 0 N − 1 x 1 ( i ) h 1 ( ( n − i ) ) N y(n)=x_1(n)\ h_1(n)=\sum_{i=0}^{N-1}{x_1(i)h_1((n-i)})_N y(n)=x1(n)⊗h1(n)=i=0∑N−1x1(i)h1((n−i))N
其中 是将 的长度扩展到N（补零）得到。
如长度为4的有限长序列x(n)={1，1，1，1},通过系统h(n)={n=0,1,2}，则有 h 1 ( n ) = 1 ， 2 ， 0 ， 0 h_1(n)={1， 2，0，0} h1(n)=1，2，0，0响应为
y ( 0 ) = x ( 0 ) h 1 ( ( 0 ) ) 4 + x 1 ( 1 ) h 1 ( ( − 1 ) ) 4 + x 1 ( 2 ) h 1 ( ( − 2 ) ) 4 + x ( 3 ) h 1 ( ( − 3 ) ) 3 y(0)=x(0)h_1((0))_4+x_1(1)h_1((-1))_4+x_1(2)h_1((-2))_4+x(3)h_1((-3))_3 y(0)=x(0)h1((0))4+x1(1)h1((−1))4+x1(2)h1((−2))4+x(3)h1((−3))3
= x 1 ( 0 ) h 1 ( 0 ) + x 1 ( 1 ) h 1 ( 3 ) + x 1 ( 2 ) h 1 ( 2 ) + x ( 3 ) h 1 ( 1 ) =x_1(0)h_1(0)+x_1(1)h_1(3)+x_1(2)h_1(2)+x(3)h_1(1) =x1(0)h1(0)+x1(1)h1(3)+x1(2)h1(2)+x(3)h1(1)
y ( 1 ) = x ( 0 ) h 1 ( ( 1 ) ) 4 + x 1 ( 1 ) h 1 ( ( 0 ) ) 4 + x 1 ( 2 ) h 1 ( ( − 1 ) ) 4 + x ( 3 ) h 1 ( ( − 2 ) ) 3 y(1)=x(0)h_1((1))_4+x_1(1)h_1((0))_4+x_1(2)h_1((-1))_4+x(3)h_1((-2))_3 y(1)=x(0)h1((1))4+x1(1)h1((0))4+x1(2)h1((−1))4+x(3)h1((−2))3
= x 1 ( 0 ) h 1 ( 1 ) + x 1 ( 1 ) h 1 ( 0 ) + x 1 ( 2 ) h 1 ( 3 ) + x ( 3 ) h 1 ( 2 ) . . . . . =x_1(0)h_1(1)+x_1(1)h_1(0)+x_1(2)h_1(3)+x(3)h_1(2)\bigm...... =x1(0)h1(1)+x1(1)h1(0)+x1(2)h1(3)+x(3)h1(2).....
参考我们上课讲到的卷积实现的算法完成下列问题
问题：序列 x ( n ) = c o s ( 0.2 π n ) x(n)=cos{(}0.2\pi n) x(n)=cos(0.2πn)，通过系统 h 1 ( n ) = 1 ， 2 ， 0 ， 0 h_1(n)={1， 2，0，0} h1(n)=1，2，0，0 ；
（1）用完成下列圆周卷积的结果（ 的长度取100）,并与的cconv(x1,x2,N)函数结果进行比较。
（2）比较圆周卷积和线性卷积（conv普通卷积）所得到的结果，分析圆周卷积的优点。（20分）
解：（1）
代码：
定义圆周卷积函数
function y=circonvt(x1,x2,N)
x_1=[x1 zeros(1,N-length(x1))];
h_1=[x2 zeros(1,N-length(x2))];
y1=conv(x_1,h_1);
z_1=[zeros(1,N) y1(1:(N-1))];
z_2=[y1((N+1):(2*N-1)) zeros(1,N)];
z=z_1(1:(2*N-1))+z_2(1:(2*N-1))+y1(1:(2*N-1));
y=z(10:N+10);
end
%% 1圆周卷积的定义如下
%（1）
clc;clear;
n=20;
t=0:99;
t1=0:n-1;
xn=cos(0.2*pi*t);
hn=[1,-2,2,-1,zeros(1,96)];
y1=circonvt(xn,hn,100);
y2=cconv(xn,hn,100);
hn=hn(1:n);
Y1=y1(1:n);
Y2=y2(1:n);
figure(1);
subplot(1,3,1);
stem(t,xn); ylabel ('hn'); xlabel ('t');title('xn原始信号'); grid on;
subplot(1,3,2);
stem(t1,Y1); ylabel ('Y1'); xlabel ('t');title('自定义圆周卷积circonvt'); grid on;
subplot(1,3,3);
stem(t1,Y2); ylabel ('Y2'); xlabel ('t');title('matlab圆周卷积cconv'); grid on;
hold on;
%%
如图1：
图1
对比上图可知，两种图像一致，说明函数构造正确。
（2）
代码：
%（2）
y3=conv(xn,hn);
Y3=y3(1:n);
figure(2);
subplot(1,3,1);
stem(t,xn); ylabel ('hn'); xlabel ('t');title('xn原始信号'); grid on;
subplot(1,3,2);
stem(t1,Y2); ylabel ('Y2'); xlabel ('t');title('matlab圆周卷积cconv'); grid on;
subplot(1,3,3);
stem(t1,Y3); ylabel ('Y3'); xlabel ('t');title('matlab线性卷积conv'); grid on;
%%
如图2：
图2
对比分析可得，圆周卷积是周期函数，线性卷积不是周期函数，线性卷积后面周期与圆周卷积一致。
2. DTFT和DFT变换
无限长序列 ，（1）画出其幅度谱；（2）将序列截取长度N=128的有限长序列，画出其DTFT频谱和DFT频谱。（3）例举计算机借助傅里叶变换分析信号频谱的步骤和方法，并指出会带来那些失真？（20分）
解：（1）
代码：
dfs函数
function[Xk]=dfs(xn,N)
n=0:N-1;
k=0:N-1;
WN=exp(-j*2*pi/N);
nk=n'*k;
Xk=(xn*WN.^nk)/N;
end
Idfs函数
function[xn]=Idfs(Xk,N)
n=0:N-1;
k=0:N-1;
WN=exp(j*2*pi/N);
nk=n'*k;
xn=(Xk*WN.^nk);
end
%% 2
% (1)
clc;clear;
n=0:0.01:2*pi;
x1n=(0.98).^n;
N = size(n,2);
Xk = dfs(x1n, N);
figure(3);
plot(n,abs(Xk));ylabel ('|Xk|'); xlabel ('n'); title('原序列的幅度谱'); %显示序列的幅度谱
hold on;
%%
如图3：
图3
（2）
代码：
%%
% (2)
n2 = 0:128;
n2 = n2*2*pi/128;
x128n=(0.98).^n2;
N2 = size(x128n,2);
X128k = dfs(x128n,N2);
k128 = 0:N2-1;
figure(4);
subplot(1,2,1);%2*2图形中的第一个
plot(k128,abs(X128k)); ylabel ('|X128k|'); xlabel ('k128'); title('128序列的DTFT频谱');
subplot(1,2,2);
stem(k128,abs(X128k)); ylabel ('|X128k|'); xlabel ('k128'); title('128序列的DFT频谱');
hold on;
%%
如图4：
图4
(3)
傅里叶变换分析信号频谱的步骤：先离散后傅里叶变换
方法：dfs /fft
由图4可知，无限长序列在截取后，其频谱会失真，此失真叫截取失真。截取长度越小失真越严重。无限长序列截短后，其频谱会产生失真，此失真叫截取失真。
3. 小波强制去噪
构造一个正弦信号 + 噪声，只采用的dwt(x,’wname’)和idwt(ca1,cd1,‘db2’)函数完成小波强制去噪的仿真。其中wname分别采用db1和db2、小波去噪采用一次分解和二次分解情况，总共4个仿真结果，分析去噪前后的信噪比。（20分）
解：
代码：
%% 3
clc;clear
X = 10*sin(0:pi/100:6*pi);
Y = awgn(X,15,'measured');
sigPower = sum(abs(X).^2)/length(X) %求出信号功率
noisePower=sum(abs(Y-X).^2)/length(Y-X) %求出噪声功率
SNR=10*log10(sigPower/noisePower) %由信噪比定义求出信噪比，单位为db
figure(5);
subplot(5,1,1);plot(Y);
xlabel('采样点'); ylabel('振幅');title(sprintf('原始信号,信噪比 %f',SNR));
% db1第一次分解
[ca1,cd1]=dwt(Y,'db1');
cd1=zeros(1,length(cd1));
X1=idwt(ca1,cd1,'db1');
X1=X1(1,1:601);
sigPower = sum(abs(X1).^2)/length(X1) %求出信号功率
noisePower=sum(abs(X1-X).^2)/length(X1-X) %求出噪声功率
SNR=10*log10(sigPower/noisePower) %由信噪比定义求出信噪比，单位为db
subplot(5,1,2);plot(X1);
xlabel('采样点'); ylabel('振幅');title(sprintf('db1第一次分解,信噪比 %f',SNR));
% db1第二次分解
[ca2,cd2]=dwt(X1,'db1');
cd2=zeros(1,length(cd2));
X2=idwt(ca2,cd2,'db1');
X2=X2(1,1:601);
sigPower = sum(abs(X2).^2)/length(X2) %求出信号功率
noisePower=sum(abs(X2-X).^2)/length(X2-X) %求出噪声功率
SNR2=10*log10(sigPower/noisePower) %由信噪比定义求出信噪比，单位为db
subplot(5,1,3);plot(X2);
xlabel('采样点'); ylabel('振幅');title(sprintf('db1第二次分解,信噪比 %f',SNR2));
% db2第一次分解
[ca3,cd3]=dwt(Y,'db2');
cd3=zeros(1,length(cd3));
X3=idwt(ca3,cd3,'db2');
X3=X3(1,1:601);
sigPower = sum(abs(X3).^2)/length(X3) %求出信号功率
noisePower=sum(abs(X3-X).^2)/length(X3-X) %求出噪声功率
SNR3=10*log10(sigPower/noisePower) %由信噪比定义求出信噪比，单位为db
subplot(5,1,4);plot(X3);
xlabel('采样点'); ylabel('振幅');title(sprintf('db2第一次分解,信噪比 %f',SNR3));
% db2第二次分解
[ca4,cd4]=dwt(X3,'db2');
cd4=zeros(1,length(cd4));
X4=idwt(ca4,cd4,'db2');
X4=X4(1,1:601);
sigPower = sum(abs(X4).^2)/length(X4) %求出信号功率
noisePower=sum(abs(X4-X).^2)/length(X4-X) %求出噪声功率
SNR4=10*log10(sigPower/noisePower) %由信噪比定义求出信噪比，单位为db
subplot(5,1,5);plot(X4);
xlabel('采样点'); ylabel('振幅');title(sprintf('db2第二次分解,信噪比 %f',SNR4));
如图5：
图5
对比前后图像可知，针对此模型型号强制去噪信噪比有了明显提高，db1比db2强制去噪明显。
4. 二次小波图像变换
用二次小波变换函数dwt2完成一幅图篇的2次小波变换，并用idwt2实现2次小波分解后的图像还原。（20分）
解：
代码：
%%
% idwt2函数
% 功能：二维离散小波反变换
clc;clear;
load woman;
nbcol = size(map,1);%返回矩阵的行数和列数
[cA1,cH1,cV1,cD1]=dwt2(X,'db1');
cod_x=wcodemat(X,nbcol);%返回矩阵X的编码矩阵，nbcol为编码的最大值
cod_cA1=wcodemat(cA1,nbcol);
cod_cH1=wcodemat(cH1,nbcol);
cod_cV1=wcodemat(cV1,nbcol);
cod_cD1=wcodemat(cD1,nbcol);
dec2d=[cod_cA1,cod_cH1;cod_cV1,cod_cD1];
A0=idwt2(cA1,cH1,cV1,cD1,'db1');
figure(10);
subplot(2,2,1),imshow(X,[]);
title('原始图像');
subplot(2,2,2),imshow(dec2d,[])
title('二维离散小波分解后的图像');
subplot(2,2,3),imshow(X,[])
title('原始图像');
subplot(2,2,4),imshow(A0,[])
title('二维小波分解重构后的图像');
%%
如图6
图6
test完整代码：
%% 1圆周卷积的定义如下
%（1）
clc;clear;
n=20;
t=0:99;
t1=0:n-1;
xn=cos(0.2*pi*t);
hn=[1,-2,2,-1,zeros(1,96)];
y1=circonvt(xn,hn,100);
y2=cconv(xn,hn,100);
hn=hn(1:n);
Y1=y1(1:n);
Y2=y2(1:n);
figure(1);
subplot(1,3,1);
stem(t,xn); ylabel ('hn'); xlabel ('t');title('xn原始信号'); grid on;
subplot(1,3,2);
stem(t1,Y1); ylabel ('Y1'); xlabel ('t');title('自定义圆周卷积circonvt'); grid on;
subplot(1,3,3);
stem(t1,Y2); ylabel ('Y2'); xlabel ('t');title('matlab圆周卷积cconv'); grid on;
hold on;
%(2)
y3=conv(xn,hn);
Y3=y3(1:n);
figure(2);
subplot(1,3,1);
stem(t,xn); ylabel ('hn'); xlabel ('t');title('xn原始信号'); grid on;
subplot(1,3,2);
stem(t1,Y2); ylabel ('Y2'); xlabel ('t');title('matlab圆周卷积cconv'); grid on;
subplot(1,3,3);
stem(t1,Y3); ylabel ('Y3'); xlabel ('t');title('matlab线性卷积conv'); grid on;%% 2
% (1)
clc;clear;
n=0:0.01:2*pi;
x1n=(0.98).^n;
N = size(n,2);
Xk = dfs(x1n, N);
figure(3);
plot(n,abs(Xk));ylabel ('|Xk|'); xlabel ('n'); title('原序列的幅度谱'); %显示序列的幅度谱
hold on;
% (2)
n2 = 0:128;
n2 = n2*2*pi/128;
x128n=(0.98).^n2;
N2 = size(x128n,2);
X128k = dfs(x128n,N2);
k128 = 0:N2-1;
figure(4);
subplot(1,2,1);%2*2图形中的第一个
plot(k128,abs(X128k)); ylabel ('|X128k|'); xlabel ('k128'); title('128序列的DTFT频谱');
subplot(1,2,2);
stem(k128,abs(X128k)); ylabel ('|X128k|'); xlabel ('k128'); title('128序列的DFT频谱');
hold on;%% 3
clc;clear
X = 10*sin(0:pi/100:6*pi);
Y = awgn(X,15,'measured');
sigPower = sum(abs(X).^2)/length(X) %求出信号功率
noisePower=sum(abs(Y-X).^2)/length(Y-X) %求出噪声功率
SNR=10*log10(sigPower/noisePower) %由信噪比定义求出信噪比，单位为db
figure(5);
subplot(5,1,1);plot(Y);
xlabel('采样点'); ylabel('振幅');title(sprintf('原始信号,信噪比 %f',SNR));
% db1第一次分解
[ca1,cd1]=dwt(Y,'db1');
cd1=zeros(1,length(cd1));
X1=idwt(ca1,cd1,'db1');
X1=X1(1,1:601);
sigPower = sum(abs(X1).^2)/length(X1) %求出信号功率
noisePower=sum(abs(X1-X).^2)/length(X1-X) %求出噪声功率
SNR=10*log10(sigPower/noisePower) %由信噪比定义求出信噪比，单位为db
subplot(5,1,2);plot(X1);
xlabel('采样点'); ylabel('振幅');title(sprintf('db1第一次分解,信噪比 %f',SNR));
% db1第二次分解
[ca2,cd2]=dwt(X1,'db1');
cd2=zeros(1,length(cd2));
X2=idwt(ca2,cd2,'db1');
X2=X2(1,1:601);
sigPower = sum(abs(X2).^2)/length(X2) %求出信号功率
noisePower=sum(abs(X2-X).^2)/length(X2-X) %求出噪声功率
SNR2=10*log10(sigPower/noisePower) %由信噪比定义求出信噪比，单位为db
subplot(5,1,3);plot(X2);
xlabel('采样点'); ylabel('振幅');title(sprintf('db1第二次分解,信噪比 %f',SNR2));
% db2第一次分解
[ca3,cd3]=dwt(Y,'db2');
cd3=zeros(1,length(cd3));
X3=idwt(ca3,cd3,'db2');
X3=X3(1,1:601);
sigPower = sum(abs(X3).^2)/length(X3) %求出信号功率
noisePower=sum(abs(X3-X).^2)/length(X3-X) %求出噪声功率
SNR3=10*log10(sigPower/noisePower) %由信噪比定义求出信噪比，单位为db
subplot(5,1,4);plot(X3);
xlabel('采样点'); ylabel('振幅');title(sprintf('db2第一次分解,信噪比 %f',SNR3));
% db2第二次分解
[ca4,cd4]=dwt(X3,'db2');
cd4=zeros(1,length(cd4));
X4=idwt(ca4,cd4,'db2');
X4=X4(1,1:601);
sigPower = sum(abs(X4).^2)/length(X4) %求出信号功率
noisePower=sum(abs(X4-X).^2)/length(X4-X) %求出噪声功率
SNR4=10*log10(sigPower/noisePower) %由信噪比定义求出信噪比，单位为db
subplot(5,1,5);plot(X4);
xlabel('采样点'); ylabel('振幅');title(sprintf('db2第二次分解,信噪比 %f',SNR4));
%% 4 二次小波变换函数dwt2
% idwt2函数
% 功能：二维离散小波反变换
clc;clear;
load woman;
nbcol = size(map,1);%返回矩阵的行数和列数
[cA1,cH1,cV1,cD1]=dwt2(X,'db1');
cod_x=wcodemat(X,nbcol);%返回矩阵X的编码矩阵，nbcol为编码的最大值
cod_cA1=wcodemat(cA1,nbcol);
cod_cH1=wcodemat(cH1,nbcol);
cod_cV1=wcodemat(cV1,nbcol);
cod_cD1=wcodemat(cD1,nbcol);
dec2d=[cod_cA1,cod_cH1;cod_cV1,cod_cD1];
A0=idwt2(cA1,cH1,cV1,cD1,'db1');
figure(6);
subplot(2,2,1),imshow(X,[]);
title('原始图像');
subplot(2,2,2),imshow(dec2d,[])
title('二维离散小波分解后的图像');
subplot(2,2,3),imshow(X,[])
title('原始图像');
subplot(2,2,4),imshow(A0,[])
title('二维小波分解重构后的图像');
5. 心得和体会
分享学习本课程的心得和体会；并与自己的研究方向相联系，所学知识如何应用到今后的科研中（1000字左右）。
现代信号处理学习心得
首先非常感谢老师这八周来的激情演讲，让我们深刻领略到现代信号处理的魅力，无论我们来自何方，专业怎么样，老师从最基础的给我们一点一点的讲解，非常难得可贵的！
回望这八周的系统性学习，深刻感受到老师的数学素养，功底十分深厚，值得我们追随学习。首先从信号的采样讲起，了解到我们这个世界基本都是把无限序列进行有限离散化处理，然后再从中发现其规律和特点，以捕捉到最细微的本质。紧接着探索卷积的本质，从手撸代码一点一点教会我们理解他是怎样的思想以及背后的数学原理。虽然这个过程中出了许多小插曲，但老师还是克服万难。像我们展示了其中的奥妙。
随机数字信号处理是由多种学科知识交叉渗透形成的，在通信、雷达、语音处理、图象处理、声学、地震学、地质勘探、气象学、遥感、生物医学工程、核工程、航天工程等领域中都离不开随机数字信号处理。随着计算机技术的进步，随机数字信号处理技术得到飞速发展。本门课主要研究了随机数字信号处理的两个主要问题:滤波器设计和频谱分析。在数字信号处理中，滤波技术占有极其重要的地位。数字滤波是语音和图像处理、模式识别、频谱分析等应用中的一个基本处理算法。但在许多应用场合，常常要处理一些无法预知的信号、噪声或时变信号，如果采用具有固定滤波系数的数字滤波器则无法实现最优滤波。在这种情况下，必须设计自适应滤波器，以使得滤波器的动态特性随着信号和噪声的变化而变化，以达到最优的滤波效果。
还学习了小波变换本文介绍了使用小波完成：间断点识别、时频分析(短时傅里叶变换)、信号去噪。下面把文中得到的一些有用结论再列举如下：小波域的高频细节的系数图像，专门用来检查原始信号间断点；随着分解级数的增加，最高级高频细节的系数图像几乎只反映间断点位置！时频分析的目的是观察原始信号中频率随时间的变换规律，窗口大小是个较为关键的参数；现实情况时频图基本上都是连在一起的；小波去噪用到的函数较多，并且很多时配套使用。
我认为接下来的学习是可以跟我的专业是相关的，基本已经确定的方向是3D打印机的故障诊断，因为这里面不仅仅还要采用视觉信息故障诊断，而且还会加入到声音震动这些性进行负荷诊断，因此小波变化是非常有用的3D打印机的故障诊断暂时分为四大类，一个是断丝，然后第二个是有层纹，第三个就是有拉丝，第四个模型就脱离平台，至于大的方面，可以采用视觉进行卷积神经网络进行区分，针对有成文和拉丝这一块儿可以采用声音传感器，进行细微的划分。再次感谢老师和同学！！！
