数学建模社区-数学中国

标题: matlab实现单边带调制及（7,4）系统循环码解码 [打印本页]

作者: 快到碗里来 时间: 2013-11-8 16:21
标题: matlab实现单边带调制及（7,4）系统循环码解码

close all; clc
clear;
echo on ;
t0=2;                %信号持续时间
ts=0.002;%抽样时间
fc=235;%载波频率
fs=1/ts;
df=0.01;%频率分辨力
t=[-t0/2:ts:t0/2];
x=cos(200*pi*t);m=x;
m(1001);
c=cos(2*pi*fc.*t);%定义载波同向分量
b=sin(2*pi*fc.*t);%定义载波正交分量
v=m.*c+imag(hilbert(m)).*b;%计算出下边带调幅分量
u=m.*c-imag(hilbert(m)).*b;%计算出上边带调幅分量
[M,m,df1]=fftseq(m,ts,df);%傅立叶变换
M=M/fs;
[U,u,df1]=fftseq(u,ts,df);%傅立叶变换
U=U/fs;
[V,v,df1]=fftseq(v,ts,df);%傅立叶变换
V=V/fs;
f=[0:df1:df1*(length(m)-1)]-fs/2;
%pause;clf;
figure(1);
subplot(2,2,1);
plot(t,m(1:length(t)));axis([-0.1,0.1,-1.1,1.1]);
xlabel('时间');title('未调信号');
%pause;
subplot(2,2,2);
plot(t,c(1:length(t)));axis([-0.1,0.1,-1.1,1.1]);
xlabel('时间');title('载波');
%pause;
subplot(2,2,3);
plot(t,u(1:length(t)));%做出上边带信号波形
axis([-0.2,0.2,-1,1.2]);
xlabel('时间');title('上边带已调信号')
%pause;
subplot(2,2,4);
plot(t,v(1:length(t)));%做出下边带波形
axis([-0.2,0.2,-1,1.2]);
xlabel('时间');title('下边带已调信号')
figure(2);
%pause;
subplot(3,1,1);
plot(f,abs(fftshift(M)));
xlabel('频率');title('模拟通信系统单边带调制，赵晓飞\2009011802/br 未调信号的频谱');
%pause;
subplot(3,1,2);
plot(f*2,abs(fftshift(U)));
title('上边带已调信号频谱');xlabel('频率');
subplot(3,1,3);plot(f,abs(fftshift(V)));
title('下边带已调信号频谱');xlabel('频率');
%%下面为对下边带信号的采样编码及解码
k=v(1:length(t));
figure(3);
subplot(2,1,1);
stem(t,k,'.');title('下边带采样信号');
axis([-0.1,0.1,-1.1,1.1]);
grid on;

n=length(k);
M=max(k);
A=(k/M)*4096;
code=zeros(n,8);

for i=1:n                         %编码
if A(i)>=0
      code(i,1)=1;
else
      code(i,1)=0;
end

if abs(A(i))>=0&&abs(A(i))<32
      code(i,2)=0;code(i,3)=0;code(i,4)=0;step=2;start=0;
elseif 32<=abs(A(i))&&abs(A(i))<64
      code(i,2)=0;code(i,3)=0;code(i,4)=1;step=2;start=32;
elseif 64<=abs(A(i))&&abs(A(i))<128
      code(i,2)=0;code(i,3)=1;code(i,4)=0;step=4;start=64;
elseif 128<=abs(A(i))&&abs(A(i))<256
      code(i,2)=0;code(i,3)=1;code(i,4)=1;step=8;start=128;
elseif 256<=abs(A(i))&&abs(A(i))<512
      code(i,2)=1;code(i,3)=0;code(i,4)=0;step=16;start=256;
elseif 512<=abs(A(i))&&abs(A(i))<1024
      code(i,2)=1;code(i,3)=0;code(i,4)=1;step=32;start=512;
elseif 1024<=abs(A(i))&&abs(A(i))<2048
      code(i,2)=1;code(i,3)=1;code(i,4)=0;step=64;start=1024;
elseif 2048<=abs(A(i))&&abs(A(i))<4096
      code(i,2)=1;code(i,3)=1;code(i,4)=1;step=128;start=2048;
end

B=floor((abs(A(i))-start)/step);
q=dec2bin(B,4)-48;
code(i,5:8)=q(1:4);
end

code=reshape(code',1,8*n);    %扩展维数
subplot(2,1,2);
stem(code,'.');axis([1 150 0 1.5]);
title('编码信号');
grid on;
%信道编码生成（7,4）系统码；
CM=xtcode(1,code);
figure(10);
subplot(2,1,1);
stem(CM,'.'); axis([1 150 0 1.5]);
title('信道编码生成系统循环码');
%将信道编码信号转换为HDB3码
CM=reshape(CM',1,length(CM)*7);
HDBM=HDB3(length(CM),CM);
subplot(2,1,2);
stem(HDBM,'-');axis([1 300 -4 4]);
title('HDB3码');
%将HDB3码解码
CDM=tranHDB3(length(HDBM),HDBM);
figure(11)
subplot(2,1,1);
stem(CDM,'-');axis([1 150 0 1.5]);
title('HDB3码解码');
%信道解码
CDM=reshape(CDM,7,length(CDM)/7);
UM=trancode(n*2,CDM');
UM=reshape(UM',1,n*8);
subplot(2,1,2);
stem(UM,'.');axis([1 150 0 1.5]);
title('信道解码输出');
%无噪声解码部分
code=reshape(code',8,n);
for p=1:n
if code(1,p)==0
      sgn(p)=-1;
else
      sgn(p)=1
end
if code(2,p)==1&&code(3,p)==1&&code(4,p)==1
      start=2048; step=128;
elseif code(2,p)==1&&code(3,p)==1&&code(4,p)==0
      start=1024;step=64;
elseif code(2,p)==1&&code(3,p)==0&&code(4,p)==1
      start=512;step=32;
elseif code(2,p)==1&&code(3,p)==0&&code(4,p)==0
      start=256;step=16;
elseif code(2,p)==0&&code(3,p)==1&&code(4,p)==1
      start=256;step=8;
elseif code(2,p)==0&&code(3,p)==1&&code(4,p)==0
      start=128;step=84;
elseif code(2,p)==0&&code(3,p)==0&&code(4,p)==1
      start=64;step=2;
elseif code(2,p)==0&&code(3,p)==0&&code(4,p)==0
      start=32;step=2;
end
hv(p)=sgn(p)*(start+step*(code(5,p)*8+code(6,p)*4+code(7,p)*2+code(8,p)*1)+step/2)/4096;
end
figure(4);
subplot(3,1,1)
stem(t,hv,'.');
title('下边带解码输出');
axis([-0.1,0.1,-1.1,1.1]);
%解调
hc1=c.*hv;
d = fdesign.lowpass(98,102,1,60,fs);
h = design(d);
hc=5*filter(h,hc1);
subplot(3,1,2);
plot(t,hc);
title('滤波后信号');
axis([-0.5,0.5,-1.1,1.1]);
[HC,hc,df1]=fftseq(hc,ts,df);
HC=HC/fs;
subplot(3,1,3);
plot(f,abs(fftshift(HC)));
title('解调输出信号频谱');
xlabel('频率');ylabel('HC')

%求傅里叶变换的子函数function [M,m,df]=fftseq(m,ts,df)fs=1/ts;if nargin==2  %nargin为输入参量的个数 n1=0;else n1=fs/df;endn2=length(m);n=2^(max(nextpow2(n1),nextpow2(n2)));%nextpow2(n)取n最接近的较大2次幂M=fft(m,n);%M为信号m的傅里叶变换，n为快速傅里叶变换的点数，及为基n-FFT变换m=[m,zeros(1,n-n2)];%构建新的m信号df=fs/n;  %重新定义频率分辨率

%先将序列变换为AMI码之后，然后再变换为HDB3码
function [HDBM]=HDB3(a,r)rn=r;l=1;                         for n=1:aif rn(n)==1 ;l=-l;                         %每检测到一个1变一次符号rn(n)=l;                      endend%AMI=rn;                            %变为AMI码：每过一个1就使原来的1变一次号m=-1;                         %符号初始化for k=1:a-3 if abs(rn(k))==1    m=rn(k);                %m记录可能存在的连续4个0前面非零码元的符号 end                               if rn(k)==0          if rn(k+1)==0             if rn(k+2)==0                   if rn(k+3)==0  %检测是否连续4个码元都是0                      rn(k+3)=2*m; %用2表示书上的V                   end             end             end       endend%rn ;                               %插入V(2)后 checked          l=1;                               %l=1表示两个V之间有偶数个非零码，0是偶数，所以初始为1，l=-1表示两个V之间有奇数个非零码for s=1:a if abs(rn(s))==2       for d=s+1:a          if abs(rn(d))==1          %V之间的非零码只有+1，-1遇到一个一使l变一次号，表示偶数或奇数             l=-l;          end          if abs(rn(d))==2             if l~=1             l=1;                      %不需要插入B时若检测到V也应把计数器清零                                           %最初版本的升级处1             else                         %检测到下一个V时，若为偶数个，则插入B             rn(d-3)=3*(-(rn(d-4))/(rn(d-4)));             %3代表B,插入不带符号的B             rn(d:a)=-rn(d:a);             end                         %V后的符号再交替          end                   end endend                               %到此处完成了插入不带符号的Bfor s=1:a if abs(rn(s))==3                               %找到B       rn(s)=rn(s)*(-1);          %B（3）与前一个位置的带符号的归一值相乘再取反，实现符号B(3)的极性与前一非符号的相反                                              %B后面第三个就是接下来的V,从它开始非零码正负号交替变化 endend                HDBM=rn;                                     %给B(3)添加了符号，并且实现了V后的符号再交替%HDB ;                                                                                  %以上便实现了HDB3的编码
%信道（7,4）循环码解码function [UM]=trancode(n,CM)UM=zeros(n,4);H=[1 0 1 1 1 0 0;1 1 1 0 0 1 0;0 1 1 1 0 0 1]; %监督矩阵SM=[0 0 0;0 0 1;0 1 0;1 0 0;1 0 1;0 1 1;1 1 1;1 1 0];%伴随式EM=[0 0 0 0 0 0 0;0 0 0 0 0 0 1;0 0 0 0 0 1 0;0 0 0 0 1 0 0;0 0 0 1 0 0 0;0 0 1 0 0 0 0;0 1 0 0 0 0 0; 1 0 0 0 0 0 0];%错误图样for k=1:n       %纠正错误 S=H'*CM(k); for p=1:7       if S==SM(p)          CM(k)=CM(k)+EM(p);                   end    end CM(k)=rem(CM(k),2); G=CM(k,:); UM(k,:)=G(:,[1,2,3,4]);end
%HDB3解码function[UM]=tranHDB3(n,HDBM)rn=HDBM;for k=1:n    if abs(rn(k))==2       rn(k-3)=0;rn(k)=0;                      %每个V（2）都变成0，V前面第三个有可能是B(3)有可能是0也恢复为0                         end if abs(rn(k))==1       rn(k)=1; endendUM=rn;

%(7,4)系统循环码，生成多项式g=（1101）%生成系统码function [CM]=xtcode(n,m)m=reshape(m',1,length(m)*n);%m=[1 0 1 1 1 0 1 1 0 0 0 0 1 1 1 1]; len=length(m); %取得m序列的长度k=rem(len,4);if k~=0 %补偿序列，使序列长度为3的倍数for p=1:3-km(len+p)=0;endend%m=reshape(m',1,8*n);m=reshape(m,4,length(m)/4);%将信息序列分组G=[1 0 0 0 1 1 0;0 1 0 0 0 1  1;0 0 1 0 1 1 1;0 0 0 1 1 0 1]; %生成矩阵CM=(m'*G);%生成系统码CM=rem(CM,2);

说明：此试验代码包括试验I和和试验II两部分，在一个试验中实现。其中ssb.m文件为主程序，其他为实现所需的函数。

欢迎光临数学建模社区-数学中国 (http://www.madio.net/)