基于MATLAB的卷积码Viterbi译码性能仿真

专业年级：电信学院09级通信与信息系统专业

小波分析课后习题

平时作业15：（选做）实现离散小波变换，绘出变换结果图。采用B2或B3滤波器，对Mallat数据进行正反离散小波变换，并绘出图形。其中：

离散小波变换的测试数据（mallat.dat）格式为ASCII码文本，内容为：

n//数据总数（整数，以下为浮点数）

x0dx//采样起点值与间隔（可用于画图）

d1//以下为n个采样数据

d2

...

Dn

%B2滤波器对256*256lena进行滤波

clc

clear all

IMG=imread('cameraman.jpg');

subplot(1,3,1);imshow(IMG);title('原图');

IMGa=double(IMG);

%%%%%%%%%%%%%%%%%%%%低通滤波%%%%%%%%%%%%%%%%%

%定义一个256*266的矩阵(加了10列0向量)

G=[zeros(256,5)IMGa zeros(256,5)];

sum=zeros(256,256);

h=[0.1250.3750.3750.125];

for k=1:length(IMG)

for j=1:length(IMG)

for i=1:length(h)

sum(k,j)=sum(k,j)+G(k,j+i-1)*h(i);

end

end

end

subplot(1,3,2);imshow(uint8(sum));title('低通滤波后图像'); %%%%%%%%%%%%%%%%%%%%高通滤波%%%%%%%%%%%%%%%%

g=[0.00610.08690.5798-0.5798-0.0869-0.0061];

clear k

clear j

clear i

SUM=zeros(256,256);

for k=1:length(IMG)

for j=1:length(IMG)

for i=1:length(g)

SUM(k,j)=SUM(k,j)+G(k,j+i-1)*g(i);

end

end

end

专业年级：电信学院09级通信与信息系统专业

subplot(1,3,3);

imshow(uint8(SUM));title('高通滤波后图像');

实验结果如下图所示:

平时作业16：（必做）实现二维哈尔小波变换的图像编码算法

程序如下:

%一维haar小波分解

function y=haardec(c,h)

%y=haardec(c,h)

y=c;

sqrt2=sqrt(2);

h=h/2;

for i=1:h

y(i)=(c(2*i-1)+c(2*i))/sqrt2;

y(h+i)=(c(2*i-1)-c(2*i))/sqrt2;

end

%二维非标准Haar小波规格化分解后图像重构程序（多级）

function a=nstdhaarrec2(c)

%输入：x载入的二维图像像数值

a=double(x);%输出：a重构后生成的图像像素数值矩阵，大小与x相同level=1;%分解重构层数

rorc=0;%作行变换（1）或列变换（2），默认为0，列行变换都做h=size(x,2);%重构的矩阵块大小，默认为整个x矩阵的变换

h1=h;

h2=h*(2^(level-1));

while h1<=h2

if rorc==1

for j=1:h1

tempcol=a(1:h1,j)';

a(1:h1,j)=haarrec(tempcol,h1)';

专业年级：电信学院09级通信与信息系统专业

end

else

if rorc==2

for i=1:h1

temprow=a(i,1:h1);

a(i,1:h1)=haarrec(temprow,h1);

end

else

for i=1:h1

temprow=a(i,1:h1);

a(i,1:h1)=haarrec(temprow,h1);

end

for j=1:h1

tempcol=a(1:h1,j)';

a(1:h1,j)=haarrec(tempcol,h1)';

end

end

end

h1+h1*2;

end

%

function y=haarrec(x,h)

%haarrec1-D haar reconstruct program

c=x;

h1=h/2;

for i=1:h1

y(2*i-1)=(c(i)+c(h1+1))/sqrt(2);

y(2*i)=(c(i)-c(h1+i))/sqrt(2);

end

%二维Haar小波非标准规格化分解程序（多级分解）

function[a,lt]=nasdhaardec2(x)

%作用:使用Haar小波对每一行和每一列像素值进行小波变换

%输入：x载入的二维图像像数值

x=imread('lena.bmp');

a=double(x);%输出：a分解后数值矩阵，大小与x相同

t=1;%记录实际分解次数

sX=size(x);

level=1;%小波分解次（级）数设定值（如果设定值超过最高可分解次数，按最高分解次数分解）

h=sX(2);%分解的矩阵块大小，默认为整个x矩阵的变换

rorc=0;%做行变换（1）或列变换（2），默认为0

lt=level;

专业年级：电信学院09级通信与信息系统专业while and(h>1,t<=level)

if rorc==1

for row=1:h

a(row,:)=haardec(a(row,:),h);

end

else

if rorc==2

for col=1:h

temp=haardec(a(:,col)',h);

a(:,col)=temp';

end

else

for row=1:h

a(row,:)=haardec(a(row,:),h);

end

for col=1:h

temp=haardec(a(:,col)',h);

a(:,col)=temp';

end

end

end

h=h/2;

t=t+1;

end

if and(h<=1,lt~=t-1)

lt=t-1;

end

codeword=codeword(simbols);

simbols=simbols-1;

simbols=simbols';

for i=1:205

Q{i,1}=simbols(i);

end

C=[Q codeword];

专业年级：电信学院09级通信与信息系统专业

平时作业17：（选做）用非标准分解方法处理256*256或更大的灰度或真彩图像（如lena），并像图10-34那样在屏幕上显示不同级别和分辨率的图像。

clc;

clear all;

IMG=imread('lena.jpg');

subplot(4,4,1);imshow(IMG);title('原图')

N=length(IMG);

dec=IMG;

dec=double(dec);

%%%%%%%%%%%%第一次行分解%%%%%%%%%%%%%%%

for m=1:N

for n=1:N/2

dec_row1(m,n)=(dec(m,2*n)+dec(m,2*n-1))/sqrt(2);dec_row1(m,n+(N/2))=(dec(m,2*n-1)-dec(m,2*n))/sqrt(2);

end

end

dec_row11=uint8(dec_row1);

subplot(4,4,2);imshow(dec_row11);title('第一次行分解')

%%%%%%%%%%%%第一次列分解%%%%%%%%%%%%%%%

for m=1:N

for n=1:N/2

dec_col1(n,m)=(dec_row1(2*n,m)+dec_row1(2*n-1,m))/sqrt(2);dec_col1(n+(N/2),m)=(dec_row1(2*n-1,m)-dec_r ow1(2*n,m))/sqrt(2);

end

end

dec_col11=uint8(dec_col1);

subplot(4,4,6);imshow(dec_col11);title('第一次列分解')

%%%%%%%%%%%%第二次行分解%%%%%%%%%%%%%%%

for m=1:N/2

N=length(IMG);

for i=1:N/4

dec_row2(m,i)=(dec_col1(m,2*i-1)+dec_col1(m,2*i))/sqrt(2);dec_row2(m,i+N/4)=(dec_col1(m,2*i-1)-dec_col1( m,2*i))/sqrt(2);

end

for i=N/2+1:length(IMG)

dec_row2(m,i)=dec_col1(m,i);

end

end

dec_row2(129:256,:)=dec_col1(129:256,:);

dec_row22=uint8(dec_row2);

专业年级：电信学院09级通信与信息系统专业

subplot(4,4,7);imshow(dec_row22);title('第二次行分解')

%%%%%%%%%%%%第二次列分解%%%%%%%%%%%%%%%%

N=length(IMG);

for m=1:N/2

N=length(IMG);

for i=1:N/4

dec_col2(i,m)=(dec_row2(2*i-1,m)+dec_row2(2*i,m))/sqrt(2);dec_col2(i+N/4,m)=(dec_row2(2*i-1,m)-dec_row2 (2*i,m))/sqrt(2);

end

for i=N/2+1:length(IMG)

dec_col2(i,m)=dec_row2(i,m);

end

end

dec_col2=[dec_col2,dec_row2(:,129:256)];

dec_col22=uint8(dec_col2);

subplot(4,4,11);imshow(dec_col22);title('第二次列分解')

%%%%%%%%%%%%第三次行分解%%%%%%%%%%%%%%%%

N=length(IMG);

for m=1:N/4

for i=1:N/8;

dec_row3(m,i)=(dec_col2(m,2*i)+dec_col2(m,2*i-1))/sqrt(2);dec_row3(m,i+N/8)=(dec_col2(m,2*i-1)-dec_col2( m,2*i))/sqrt(2);

end

for i=N/4+1:length(IMG)

dec_row3(m,i)=dec_col2(m,i);

end

end

dec_row3=[dec_row3;dec_col2(65:256,:)];

dec_row33=uint8(dec_row3);

subplot(4,4,12);imshow(dec_row33);title('第三次行分解')

%%%%%%%%%%%%第三次列分解%%%%%%%%%%%%%%%%%

N=length(IMG);

for m=1:N/4

for i=1:N/8;

dec_col3(i,m)=(dec_row3(2*i,m)+dec_row3(2*i-1,m))/sqrt(2);dec_col3(i+N/8,m)=(dec_row3(2*i-1,m)-dec_row3 (2*i,m))/sqrt(2);

end

for i=N/4+1:length(IMG)

dec_col3(i,m)=dec_row3(i,m);

end

专业年级：电信学院09级通信与信息系统专业

end

dec_col3=[dec_col3,dec_row3(:,65:256)];subplot(4,4,16)

dec_col33=uint8(dec_col3);

subplot(4,4,16);imshow(dec_col33);title('第三次列分解

')

平时作业18：（选做）用整数小波的S 或2/6变换对256*256Lena 灰度图像进行非标准方法的3级分解与重构。

1.总体设计

用下面的已知分解公式进行小波分解和重构。分解公式：

1,,21,21,,21,[]

j k j k j k j k j k j k d s s s s d -+--=-=+重构公式：

,21,1,,211,,2[]j k j k j k j k j k j k

s s d s d s --+-=-=+，其中[]表示取整。

进行非标准小波分解，即交替进行3次行变换和3次列变换，程序对每次变换后的结果都保存为位图文件，运行后可以在程序所在路径下看到保存的6个分解位图文件和6个重构位图文件。2.实现方法

编写S 变换的分解和重构子程序，分别对图像数据进行一次行列分解和列行重构，程序返回该次变换后的行列矩阵，在主程序中可以连续三次调用行列变换，即完成对原始图像的3级分解和重构，这里的变换是完全可逆的，也就是能够完全恢复原图像数据。程序代码如下:clc;clear all;

IMG=imread('lena.jpg');

subplot(2,3,1);imshow(IMG);title('原图')N=length(IMG);

专业年级：电信学院09级通信与信息系统专业

R0=double(IMG); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%第一级的奇数序列%%%%%%%%%%%%%%%%

for m=1:N

for i=1:2:N-1

O1(m,((i+1)/2))=R0(m,i);

end

end

%得到第一级偶序列

for m=1:N

for i=2:2:N

E1(m,i/2)=R0(m,i);

end

end

i=0;

for m=1:N

for i=1:N/2

%各行差值数据；不对细节数据再做变换，只对S1再预测

d1(m,i)=(O1(m,i)-E1(m,i))/2;

%更新过程

S1(m,i)=E1(m,i)+d1(m,i)/2;

end

end

L1=[S1,d1];

L18=uint8(L1);

M=length(IMG);

%得到列变换的第一级奇序列

for m=1:N

for i=1:2:N-1

O11(((i+1)/2),m)=L1(i,m);

end

end

%得到列变换的第一级偶序列

for m=1:N

for i=2:2:N

E11(i/2,m)=L1(i,m);

end

end

i=0;

for m=1:N

for i=1:N/2

%各行差值数据；不对细节数据再做变换，只对S1再预测

d11(i,m)=(O11(i,m)-E11(i,m))/2;

专业年级：电信学院09级通信与信息系统专业

%更新过程

S11(i,m)=E11(i,m)+d11(i,m)/2;

end

end

L11=[S11;d11];

L118=uint8(L11);

subplot(2,3,2);imshow(L118);title('整数小波第一级分解') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%第二级的奇数序列%%%%%%%%%%%%%%%%

B=L11(1:128,1:128);

i=0;

%得到第二级奇序列

for m=1:N/2

for i=1:2:N/2-1

O2(m,((i+1)/2))=B(m,i);

end

end

i=0;

%得到第二级偶序列

for m=1:N/2

for i=2:2:N/2

E2(m,i/2)=B(m,i);

end

end

i=0;

for m=1:N/2

for i=1:N/4

%各行差值数据；不对细节数据再做变换，只对S1再预测

d2(m,i)=(O2(m,i)-E2(m,i))/2;

%更新过程

S2(m,i)=E2(m,i)+d2(m,i)/2;

end

end

L2=[S2,d2];%完成了第二次行分解

%第二级列的奇数序列

i=0;

for m=1:N/2

for i=1:2:N/2-1

O22(((i+1)/2),m)=L2(i,m);

end

end

i=0;

%得到第二级列的偶序列

专业年级：电信学院09级通信与信息系统专业

for m=1:N/2

for i=2:2:N/2

E22(i/2,m)=L2(i,m);

end

end

i=0;

for m=1:N/2

for i=1:N/4

%各列差值数据；不对细节数据再做变换，只对S11再预测

d22(i,m)=(O22(i,m)-E22(i,m))/2;

%更新过程

S22(i,m)=E22(i,m)+d22(i,m)/2;

end

end

C=[S22;d22];Y=[C,L11(1:128,129:256)];L22=[Y;d11];

subplot(2,3,3);imshow(uint8(L22));title('整数小波第二级分解') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%第三级的奇数序列%%%%%%%%%%%%%%%%

D=L22(1:64,1:64);

i=0;

%得到第三级奇序列

for m=1:N/4

for i=1:2:N/4-1

O3(m,((i+1)/2))=D(m,i);

end

end

i=0;

%得到第三级偶序列

for m=1:N/4

for i=2:2:N/4

E3(m,i/2)=D(m,i);

end

end

i=0;

for m=1:N/4

for i=1:N/8

%各行差值数据；不对细节数据再做变换，只对S1再预测

d3(m,i)=(O3(m,i)-E3(m,i))/2;

%更新过程

S3(m,i)=E3(m,i)+d3(m,i)/2;

end

end

L3=[S3,d3];

专业年级：电信学院09级通信与信息系统专业

%得到列变换的第一级奇序列

for m=1:N/4

for i=1:2:N/4-1

O33(((i+1)/2),m)=L3(i,m);

end

end

%得到列变换的第一级偶序列

for m=1:N/4

for i=2:2:N/4

E33(i/2,m)=L3(i,m);

end

end

i=0;

for m=1:N/4

for i=1:N/8

%各行差值数据；不对细节数据再做变换，只对S1再预测

d33(i,m)=(O33(i,m)-E33(i,m))/2;

%更新过程

S33(i,m)=E33(i,m)+d33(i,m)/2;

end

end

x=[S33;d33];x=[x,L22(1:64,65:128)];x=[x;d22];

y=[x,L11(1:128,129:256)];z=[y;d11];

subplot(2,3,4);imshow(uint8(z));title('整数小波第三级分解') subplot(2,3,5);imshow(L118(1:128,1:128));title('整数小波重构')

专业年级：电信学院09级通信与信息系统专业

平时作业19：（选做）用Le Gall5/3整数小波或Daubechies9/7离散小波的提升实现，对256*256Lena灰度图像进行非标准方法的3级分解与重构。

%=====该程序对向量进行小波变换，所用的小波为利用提升方案实现的Daubechies9/7小波===== function[cn,dn]=lifting_db97(x)

%这是一维小波分解

%cn和dn分别是分解得到的低频向量(概貌系数)和高频向量(细节系数)

%x为待分解的向量

L=floor(length(x)/2);

s0=zeros(1,L);d0=zeros(1,L);

alpha=-1.586134342;beta=-0.05298011854;gamma=0.8829110762;delta=0.4435068522;

K=1.149604398;%设定参数值

for i=1:L%进行奇偶二抽取

s0(i)=x(2*i-1);

d0(i)=x(2*i);

end

s1=s0(2:L);

s1=[s1,s0(L)];

d0=d0+alpha*(s0+s1);

d1=d0(1:L-1);

d1=[d0(1),d1];

s0=s0+beta*(d0+d1);

s1=s0(2:L);

s1=[s1,s0(L)];

d0=d0+gamma*(s0+s1);

d1=d0(1:L-1);

d1=[d0(1),d1];

s0=s0+delta*(d0+d1);

cn=K*s0;dn=d0/K;%求出cn和dn

if L

dn=[dn,dn(L)];

end

对二维信号（例如二维图像），实现多级小波的分解和重构，用到的小波函数是DBN小波,即Daubechies 小波

function X=fn_DBn_2D(X,wname,level)

%size(X):(2^n)*(2^n)

%wname:'dbN',N=1,2,3,4,5,...

[L,H,L_r,H_r]=wfilters(wname);

专业年级：电信学院09级通信与信息系统专业lf=length(L);

if size(X,1)~=size(X,2),error('square image needed');end

N=size(X,1);

if level>0

for ii=1:level%%%%%%%%%%%%%%%%%行变换

half=N/2;

for m=1:N

i=1;j=1;

while j<=N-(lf-1)

temp(m,i)=0;

temp(m,i+half)=0;

for k=1:lf

jk=j+k-1;

temp(m,i)=temp(m,i)+X(m,jk)*L(k);

temp(m,i+half)=temp(m,i+half)+X(m,jk)*H(k);

end

i=i+1;j=j+2;

end

while j<=N-1

temp(m,i)=0;

temp(m,i+half)=0;

for k=1:lf

jk=j+k-1;

if jk>N,jk=jk-N;end

temp(m,i)=temp(m,i)+X(m,jk)*L(k);

temp(m,i+half)=temp(m,i+half)+X(m,jk)*H(k);

end

i=i+1;j=j+2;

end

for i=1:N

X(m,i)=temp(m,i);

end

end

for m=1:N%%%%%%%%%%%%%%%%%列变换i=1;j=1;

while j<=N-(lf-1)

temp(i,m)=0;

temp(i+half,m)=0;

for k=1:lf

jk=j+k-1;

temp(i,m)=temp(i,m)+X(jk,m)*L(k);

temp(i+half,m)=temp(i+half,m)+X(jk,m)*H(k);

end

专业年级：电信学院09级通信与信息系统专业

i=i+1;j=j+2;

end

while j<=N-1

temp(i,m)=0;

temp(i+half,m)=0;

for k=1:lf

jk=j+k-1;

if jk>N,jk=jk-N;end

temp(i,m)=temp(i,m)+X(jk,m)*L(k);

temp(i+half,m)=temp(i+half,m)+X(jk,m)*H(k);

end

i=i+1;j=j+2;

end

for i=1:N

X(i,m)=temp(i,m);

end

end

N=N/2;

end

else

level=-level;

half=N/2^level;

temp=zeros(N,N);

for jj=1:level

for m=1:2*half

j=1;i=half-floor((lf-1)/2)+1;

while j<(lf-1)

temp(j,m)=0;

temp(j+1,m)=0;

for k=1:2:(lf-1)

ik=i+floor(k/2);

if ik>half,ik=ik-half;end

temp(j,m)=temp(j,m)+X(ik,m)*L_r(k+1)+X(ik+half,m)*H_r(k+1);

temp(j+1,m)=temp(j+1,m)+X(ik,m)*L_r(k)+X(ik+half,m)*H_r(k);

end

j=j+2;

i=i+1;

end

j=lf-1;i=1;

while(i+lf/2-1)<=half

temp(j,m)=0;

temp(j+1,m)=0;

for k=1:2:(lf-1)

专业年级：电信学院09级通信与信息系统专业

ik=i+floor(k/2);

temp(j,m)=temp(j,m)+X(ik,m)*L_r(k+1)+X(ik+half,m)*H_r(k+1);

temp(j+1,m)=temp(j+1,m)+X(ik,m)*L_r(k)+X(ik+half,m)*H_r(k);

end

j=j+2;

i=i+1;

end;

for i=1:2*half

X(i,m)=temp(i,m);

end

end

for m=1:2*half

j=1;i=half-floor((lf-1)/2)+1;

while j<(lf-1)

temp(m,j)=0;

temp(m,j+1)=0;

for k=1:2:(lf-1)

ik=i+floor(k/2);

if ik>half,ik=ik-half;end

temp(m,j)=temp(m,j)+X(m,ik)*L_r(k+1)+X(m,ik+half)*H_r(k+1);

temp(m,j+1)=temp(m,j+1)+X(m,ik)*L_r(k)+X(m,ik+half)*H_r(k);

end

j=j+2;

i=i+1;

end

j=lf-1;i=1;

while(i+lf/2-1)<=half

temp(m,j)=0;

temp(m,j+1)=0;

for k=1:2:(lf-1)

ik=i+floor(k/2);

temp(m,j)=temp(m,j)+X(m,ik)*L_r(k+1)+X(m,ik+half)*H_r(k+1);

temp(m,j+1)=temp(m,j+1)+X(m,ik)*L_r(k)+X(m,ik+half)*H_r(k);

end

j=j+2;

i=i+1;

end;

for i=1:2*half

X(m,i)=temp(m,i);

end

end

half=half*2;

end

专业年级：电信学院09级通信与信息系统专业

end

平时作业20：实现EZW,SPIHT算法。

EZW算法原理

给定一幅图像,图像压缩程序将对他进行二维离散小波变换。将变换结果传给EZW算法,量化器对其进行量化,根据公式T0=2[log2(k)]求得初始阈值其中T0初始阈值,k为系数的最大绝对值。由此可以看出,初始阈值为小于系数最大绝对值2的最大整数次方。运用初始阈值,对所有系数进行主扫描,确定哪些是绝对值大于初始阈值即重要的。完成主扫描后,当前阈值减半,然后进行辅扫描,它是对重要系数的细节信息进行渐进编码。这个过程(从主扫描开始)反复进行直到希望的精度实现为止。

希望的精度是指阈值达到最小值或者是预定的性能指标实现。最常用的性能指标是压缩后图像的每像素比特数,当达到预定的每像素比特数就停止编码。因为每像素比特数是可以可控的,这可以确保准确的压缩比。另一个性能指标是编码误差,追踪原始图像和重构图像间最大绝对值像素误差,当误差落入预定范围时停止编码。