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文件为主程序，其他为实现所需的函数。