《基于QAM调制的无线衰落信道的性能分析与仿真毕业设计(32页).doc》由会员分享，可在线阅读，更多相关《基于QAM调制的无线衰落信道的性能分析与仿真毕业设计(32页).doc（32页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、-基于QAM调制的无线衰落信道的性能分析与仿真毕业设计-第 30 页 实践教学 兰州理工大学计算机与通信学院2013年秋季学期通信系统综合训练 题 目：基于QAM调制的无线衰落信道的性能分析与仿真 专业班级： 通信工程（1）班 姓 名： 学 号： 指导教师： 成 绩： 摘要近年来，信息通信领域中，发展最快、应用最广的就是无线通信技术。但无线信道中的衰落现象，严重影响通信系统的性能，所以了解和掌握衰落信道中无线通信系统的性能是重中之重。本次课题正是基于QAM调制的无线衰落信道的性能分析与仿真，首先介绍了QAM调制解调原理以及无线衰落信道特征，其次利用Matble分析工具分析16QAM已调信号的频
2、谱图，最后分析衰落信道下的系统的误码率，并与高斯信道下的性能进行对比。关键词：正交振幅调制；Matble；调制解调；仿真目录摘要1前言3一 基本原理41.1QAM调制解调原理41.1.1QAM调制原理41.1.2QAM的解调和判决原理51.1.3QAM的误码性能6二 无线衰落信道的特征82.1瑞利衰落信道82.2瑞利衰落信道基本模型8三 系统分析10四 系统设计及系统调试114.1 16QAM调制信号114.1.1 信号源114.1.2 串并转换114.1.3 2-4电平转换124.1.4 增加载波134.1.5 调制信号形成144.2 16QAM调制信号的噪声叠加154.3 16QAM解调模
3、块的建立与仿真154.3.1 滤波器164.3.2 抽样判决和4-2电平转换164.3.3 并串转换174.4 无线衰落信道性能分析184.4.1 16QAM抗噪声性能仿真18总结20参考文献21附录22致谢36前言随着通信业迅速的发展，传统通信系统的容量已经越来越不能满足当前用户的要求，而可用频谱资源有限，业不能靠无限增加频道数目来解决系统容量问题。另外，人们亦不能满足通信单一的语音服务，希望能利用移动电话进行图像等多媒体信息的通信。但由于图像通信比电话需要更大的信道容量。高效、可靠的数字传输系统对于数字图像通信系统的实现很重要，正交幅度调制QAM是数字通信中一种经常利用的数字调制技术，尤其
4、是多进制QAM具有很高的频带利用率，在通信业务日益增多使得频带利用率成为主要矛盾的情况下，正交幅度调制方式是一种比较好的选择。如今每一天大约有15万人成为新的无线用户，全球范围内的无线用户数量目前已经超过2亿。这些人包括大学教授、仓库管理员、护士、商店负责人、办公室经理和卡车司机。他们使用无线通信技术的方式和他们自身的工作一样都在不断地更新。但是在无线信道中存在着衰落现象，这将严重影响通信系统的性能。所以了解和掌握衰落信道中无线通信系统的性能成为一个关键问题。一 基本原理1.1 QAM调制解调原理1.1.1 QAM调制原理正交幅度调制（QAM）是数字通信中一种经常利用的数字调制技术，尤其是多进
5、制QAM具有很高的频带利用率，在通信业务日益增多使得频带利用率成为主要矛盾的情况下，正交幅度调制方式是一种比较好的选择。正交幅度调制（QAM）信号采用了两个正交载波，每一个载波都被一个独立的信息比特序列所调制。发送信号波形如图1-1所示： （1-1）图1-1 M=16QAM信号星座图式（1-1）中和是电平集合，这些电平是通过将k比特序列映射为信号振幅而获得的。例如一个16位正交幅度调制信号的星座图如下图所示，该星座是通过用M4PAM信号对每个正交载波进行振幅调制得到的。利用PAM分别调制两个正交载波可得到矩形信号星座。QAM 可以看成是振幅调制和相位调制的结合。因此发送的QAM信号波形可表示为
6、（1-2）式 （1-2）如果那么QAM方法就可以达到以符号速率同时发送个二进制数据。图1-2给出了QAM调制器的框图。图1-2 QAM调制器框图1.1.2 QAM的解调和判决原理假设在信号传输中存在载波相位偏移和加性高斯噪声。因此r(t)可以表示（1-3）式： （1-3）其中是载波相位偏移，且n(t)表示为（1-4）所示： （1-4）将接收信号与下述两个相移函数式（1-5）和（1-6）进行相关运算 （1-5） （1-6）如图1-3所示，相关器的输出抽样后输入判决器。使用图1-3中所示的锁相环估算接收信号的载波相位偏移，相移和对该相位偏移进行补偿。图1-3 QAM信号的解调和判决假设图中所示的时
7、钟与接收信号同步，以使相关器的输出在适当的时刻及时被抽样。在这些条件下两个相关器的输出分别表示为式（1-7）和式（1-8）所示： （1-7） （1-8）其中nc和ns分别表示为式（1-9）和式（1-10）所示： （1-9） （1-10）噪声分量是均值为0，方差为的互不相关的高斯随机变量。最佳判决器计算距离量度表示为（1-11）式所示： （1-11）1.1.3 QAM的误码性能矩形QAM信号星座最突出的优点就是容易产生PAM信号可直接加到两个正交载波相位上，此外它们还便于解调。对于M下的矩形信号星座图（k为偶数），QAM信号星座图与正交载波上的两个PAM信号是等价的，这两个信号中的每一个上都有个
8、信号点。因为相位正交分量上的信号能被相干判决极好的分离，所以易于通过PAM的误码率确定QAM的误码率。M进制QAM系统正确判决的概率表示为（1-12）所示： （1-12）其中（1-12）式中是进制PAM系统的误码率，该PAM系统具有等价QAM系统的每一个正交信号中的一半平均功率。通过适当调整M进制PAM系统的误码率，可得式（1-13）所示： （1-13）其中是每个符号的平均信噪比。因此M进制QAM的误码率为（1-14）所示： ) （1-14）可以注意到，当k为偶数时，这个结果对M情形时精确的，而当k为奇数时，就找不到等价的进制PAM系统。如果使用最佳距离量度进行判决的最佳判决器，可以求出任意k
9、1误码率的严格上限表示为（1-15）所示： （1-16）其中是每比特的平均信噪比。二 无线衰落信道的特征2.1瑞利衰落信道在陆地移动通信中，移动台往往受到各种障碍物和其他移动体的影响，以致到达移动台的信号是来自不同传播路径的信号之和。而描述这样一种信道的常用信道模型便是瑞利衰落信道。瑞利衰落信道（Rayleigh fading channel）是一种无线电信号传播环境的统计模型。这种模型假设信号通过无线信道之后，其信号幅度是随机的，表现为“衰落”特性，并且多径衰落的信号包络服从瑞利分布。由此，这种多径衰落也称为瑞利衰落。 这一信道模型能够描述由电离层和对流层反射的短波信道，以及建筑物密集的城市
10、环境。瑞利衰落只适用于从发射机到接收机不存在直射信号的情况，否则应使用莱斯衰落信道作为信道模型。假设经反射（或散射）到达接收天线的信号为N个幅值和相位均随机的且统计独立的信号之和。信号振幅为r,相位为,则其包络概率密度函数为（2-1）所示：P(r)= (r0) （2-1）相位概率密度函数为（2-2）所示：P()=1/2 （） （2-2）2.2瑞利衰落信道基本模型根据ITU-RM.1125标准，离散多径衰落信道模型如下，可表示为式（2-3）所示： （2-3）其中复路径衰落，服从瑞利分布；是多径时延。多径衰落信道模型框图如图2-1所示：图2-1 多径衰落信道模型框图假设经反射（或散射）到达接收天线
11、的信号为N个幅值和相位均随机的且统计独立的信号之和。信号振幅为r,相位为,则信号经过瑞利衰落信道其包络概率密度函数为（2-4）所示：P(r)= (r0) （2-4）相位概率密度函数为（2-5）所示：P()=1/2（） （2-5）三 系统分析本次课题是基于QAM调制的无线衰落信道的性能分析与仿真，首先对信号源进行16QAM调制，并分析已调信号的频谱；其次让信号通过瑞利信道并进行QAM解调；最后分析了衰落信道下的系统的误码率，其系统模型如下图3-1所示：信号源QAM调制瑞利衰落信道受信者误码性能分析已调信号频谱QAM解调噪声源图3-1 系统分析模型四 系统设计及系统调试4.1 16QAM调制信号4
12、.1.1 信号源本程序中，信号源为8位二进制代码-1 1 -1 1 1 1 -1 -1， sigexpand函数的作用是将代码扩展为码元宽度为1的双极性波形，如下图4-1所示： 图4-1 二进制代码及星座图4.1.2 串并转换信号源通过串并变换，将原来的一路信源信号变成两路信号，分别为上支路信号和下支路信号，独立地进行调制和解调。串并变换的规则是根据序列编号的奇偶行，将编号为奇的码元编成一路信号，将编号为偶的码元编成一路信号。经过串并转换后，并行输出的每一路码元传输速率降为原来的一半即Rb/2. 输入d：-1 1 -1 1 1 1 -1 -1 上支路d_NRZ1：-1 -1 1 -1 下支路d
13、_NRZ2：1 1 1 -1图4-2 串并转换后上下支路信号时域波形图4.1.3 2-4电平转换2-4电平转换就是将输入信号的2电平信号状态经过转换后变成相应的4电平信号。这里选择的映射关系如下所示：映射前数据 双极性 电平/V 00 -1 -1 -3 01 -1 1 -1 10 1 -1 1 11 1 1 3根据以上的映射关系，可得到上下支路分别为上支路d_NRZ1：-1 -1 1 -1；下支路d_NRZ2：1 1 1 -1；2-4电平转换信号：上支路d_NRZ1：-1 -1 1 -12-4电平转换后：-3 1下支路d_NRZ2：1 1 1 -12-4电平转换后：3 1 图4-3 2-4电平
14、转换后上下支路信号时域波形图这里4电平信号的码元传输速率已降为Rb/44.1.4 增加载波在本课题中，选用的载波是载波幅度A=1，载波频率fc=2Hz,上支路分量的载波是h1t=A*cos(2*pi*fc*t),正交分量的载波是h2t=A*sin(2*pi*fc*t) 。上下支路信号在加载波之前还经过平滑处理，以滤除较高频率的信号，使实验结果更加理想。上下支路信号加载波后的图形如图4-4所示：图4-4上下支路调制信号时域波形图4.1.5 调制信号形成上下支路调制信号形成后，将两个分量相加，既可得到16QAM调制信号，如下图4-5所示：图4-5 已调信号波形图4.2 16QAM调制信号的噪声叠加
15、 本次仿真采用的噪声是高斯白噪声，这是一种最常见的噪声，白噪声的功率谱密度在所有频率上均为一常数，且仅在t=0时才相关，而在任意两个时刻的随机变量都是不相关的。 对已调制信号可采用wgn函数添加加性高斯噪声。y = wgn(m,n,p) 产生一个m行n列的高斯白噪声的矩阵，p以dBW为单位指定输出噪声的强度。为使解调效果较好，采用噪声的强度较小，设置Pn=-10dB.4.3 16QAM解调模块的建立与仿真 系统先前所得的16QAM调制信号通过高斯白噪声信道以后便可以解调了。本文所采用的解调器原理为相干解调法，即已调信号与载波相乘，送入到低通滤波器，其对应原理图中信号输入并与载波相乘后通过LPF
16、的部分，输出送入到判决器判决，再经4-2电平转换和并串转换即可得到解调信号。4.3.1 滤波器 IIR滤波器采用的巴特沃斯低通滤波器有现成的模型，我们可以加以利用，因此在本文涉及的仿真中滤波器均选择巴特沃斯低通滤波器。巴特沃斯低通滤波器的标准形式为：b,a=butter(N,W0); y=filter(b,a,x); N表示要选取的低通滤波器的阶数，W0表示滤波器的截止频率，b,a为滤波器返回的特性参数。在第二行程序中，x表示输入序列，y表示输出序列。整个函数表示信号通过滤波器的过程。 图4-6上下支路通过低通滤波器信号时域波形图4.3.2 抽样判决和4-2电平转换 抽样判决是在每个码元中间抽
17、样，并用几个判决语句进行判决，例如当码元幅度大于1时，判为3。4-2电平转换是2-4电平转换的逆过程，其映射关系如下图4-7所示：映射前数据 电平/V 双极性 -3 00 -1 -1 -1 01 -1 1 1 10 1 -1 3 11 1 1 图4-7上下支路抽样判决及4-2转换后信号时域波形图4.3.3 并串转换 经过并串转换即可得到解调信号，采用循环语句实现，for s=1:N/2 ddd(2*s-1)=dd111111(s); ddd(2*s)=dd222222(s);end将两路并联信号经过转换成为一路信号，下图为基带信号和解调信号的对比及它们的频率谱密度图如图4-8所示：图4-8基带
18、与解调信号可以看出，并没有错码，可见本次仿真是成功的。4.4 无线衰落信道性能分析4.4.1 16QAM抗噪声性能仿真对于QAM，可以看成是由两个相互正交且独立的多电平ASK信号叠加而成。因此，利用多电平误码率的分析方法，可得到M进制QAM的误码率为： (4-1-1) 式中，Eb为每码元能量，n0为噪声单边功率谱密度。通过调整高斯白噪声信道的信噪比snr（Eb/No），可以得到如图4-9所示的误码率图： 图4-9 QAM信号误码率分析可见16QAM信号的误码率随着信噪比的增大而逐渐减小，这与理论分析是完全一致的。 总结通信综合训练是培养学生综合运用所学的理论知识，发现、提出、分析和解决实际问题
19、，锻炼实践能力的重要环节，是对学生实际工作能力的具体训练和考察过程。本次课题是基于QAM调制的无线衰落信道的性能分析与仿真，首先介绍了QAM调制解调原理以及无线衰落信道特征，其次利用Matble分析工具分析16QAM已调信号的频谱图，最后分析衰落信道下的系统的误码率，并与高斯信道下的性能进行对比。在设计过程中困难有很多，其主要表现在不熟练软件编程。总之，通过此次设计，我的收获有以下几点：首先，通过此次课程设计我掌握了无线衰落信道的传输特征，并且在实际操作中掌握了QAM的调制解调原理；其次，通过此次仿真训练我清楚地了解了16QAM已调信号的频谱特征；最后，我觉得细节决定成败。不做系统，许多细小的
20、环节是注意不到的，而这诸多环节往往影响你整个系统的正常运转。这可真应验了“细节决定一切”这句话。这一切告诉我做任何事情必须从全局出发，并且要注意其中的任何一个细节。参考文献1樊昌信，曹丽娜编著通信原理北京：国防工业出版社，20102樊昌信通信原理北京：国防工业出版社，20023曹志刚等著现代通信原理北京：清华大学出版社，20014吴伟陵等著移动通信原理北京：电子工业出版社，20055.李建新现代通信系统分析与仿真-MATLAB 通信工具箱西安：西安电子科技大学出版社，20006潘子宇Matlab通信仿真设计指导书南京工程学院，20117刘敏MATLAB通信仿真与应用北京：国防工业出版社附录1、
21、将输入的序列扩展成间隔为N-1个0的序列function out=sigexpand(d,M)N=length(d);out=zeros(M,N);out(1,:)=d;out=reshape(out,1,M*N);2、计算信号的傅里叶变化functionf,sf=T2F(t,st);dt=t(2)-t(1);T=t(end);df=1/T;N=length(st);f=-N/2*df:df:N/2*df-df;sf=fft(st);sf=T/N*fftshift(sf);4.1.2 主函数代码fc=2; % 载波频率N_sample=8; % 基带码元抽样点数N=8; % 码元数Ts=1;
22、% 码元宽度dt=Ts/fc/N_sample; % 抽样时间间隔T=N*Ts; % 信号持续时间长度t=0:dt:T-dt; % 时间向量Lt=length(t); % 时间向量长度tx1=0; % 时域波形图横坐标起点tx2=8; %时域波形图横坐标终点ty1=-4.5; %时域波形图纵坐标起点ty2=4.5; %时域波形图纵坐标终点fx1=-10; %功率谱图横坐标起点fx2=10; %功率谱图横坐标终点fy1=-40; %功率谱图纵坐标起点fy2=25; %功率谱图纵坐标终点%产生二进制信源d=-1,1,-1,1,1,1,-1,-1dd=sigexpand(d,fc*N_sample)
23、; % 双极性gt=ones(1,fc*N_sample); % NRZ波形d_NRZ=conv(dd,gt); % 基带信号 figure(1);subplot(2,3,1);plot(t,d_NRZ(1:Lt);axis(tx1,tx2,ty1,ty2);xlabel(时间(S);ylabel(幅度);title(基带信号时域波形图);grid;figure(5);subplot(2,2,1);plot(t,d_NRZ(1:Lt);axis(tx1,tx2,ty1,ty2);xlabel(时间(S);ylabel(幅度);title(基带信号时域波形图);grid;f1,d_NRZf=T2
24、F(t,d_NRZ(1:Lt);figure(5);subplot(2,2,2);plot(f1,10*log10(abs(d_NRZf).2/T);axis(fx1,fx2,fy1,fy2);xlabel(频率(Hz);ylabel(功率谱密度(dB/Hz);title(基带信号功率谱图);grid;% 串并转换d1=;d2=;for i=1:N/2 d1(i)=d(2*(i-1)+1); d2(i)=d(2*(i-1)+2);end% 上支路dd1=sigexpand(d1,2*fc*N_sample); gt1=ones(1,2*fc*N_sample); d_NRZ1=conv(dd1
25、,gt1); figure(1);subplot(2,3,2);plot(t,d_NRZ1(1:Lt);axis(tx1,tx2,ty1,ty2);xlabel(时间(S);ylabel(幅度);title(串并转换后上支路信号时域波形图);grid;% 下支路dd2=sigexpand(d2,2*fc*N_sample); d_NRZ2=conv(dd2,gt1); figure(1);subplot(2,3,5);plot(t,d_NRZ2(1:Lt);axis(tx1,tx2,ty1,ty2);xlabel(时间(S);ylabel(幅度);title(串并转换后下支路信号时域波形图);
26、grid;% 载波h1t=cos(2*pi*fc*t);h2t=sin(2*pi*fc*t);figure(1);subplot(2,3,4);plot(t,h1t);axis(tx1,tx2,ty1,ty2);xlabel(时间(S);ylabel(幅度);title(载波信号时域波形图);grid;% 上下支路2-4电平转换d11=; % 上支路d22=; % 下支路for m=1:N/4; d11(m)=2*d1(2*m-1)+d1(2*m); d22(m)=2*d2(2*m-1)+d2(2*m);enddd11=sigexpand(d11,4*fc*N_sample); % 上支路gt
27、2=ones(1,4*fc*N_sample); d_NRZ11=conv(dd11,gt2); figure(1);subplot(2,3,3);plot(t,d_NRZ11(1:Lt);axis(tx1,tx2,ty1,ty2);xlabel(时间(S);ylabel(幅度);title(2-4电平转换后上支路信号时域波形图);grid;figure(3);subplot(2,2,1);plot(t,d_NRZ11(1:Lt);axis(tx1,tx2,ty1,ty2);xlabel(时间(S);ylabel(幅度);title(2-4转换后上支路信号时域波形图);grid;dd22=si
28、gexpand(d22,4*fc*N_sample); % 下支路d_NRZ22=conv(dd22,gt2); figure(1);subplot(2,3,6);plot(t,d_NRZ22(1:Lt);axis(tx1,tx2,ty1,ty2);xlabel(时间(S);ylabel(幅度);title(2-4转换后下支路信号时域波形图);grid;figure(3);subplot(2,2,3);plot(t,d_NRZ22(1:Lt);axis(tx1,tx2,ty1,ty2);xlabel(时间(S);ylabel(幅度);title(2-4转换后下支路信号时域波形图);grid;%
29、 对上下支路信号进行平滑pd_NRZ11=d_NRZ11(1:Lt);pd_NRZ22=d_NRZ22(1:Lt);tao=3/16;tm=0:dt:8*Ts*tao-dt;for w=2:N/4 phxh1=d11(w-1)+(d11(w)-d11(w-1)*0.5*(1+sin(pi*tm/8*tao*Ts); phxh2=d22(w-1)+(d22(w)-d22(w-1)*0.5*(1+sin(pi*tm/8*tao*Ts); for k=1:8*tao*fc*N_sample pd_NRZ11(w-1)*64+k)=phxh1(k); pd_NRZ22(w-1)*64+k)=phxh2
30、(k); endend% 生成16QAM信号s_16qam1=pd_NRZ11 .* h1t; %生成上支路频带信号s_16qam2=pd_NRZ22 .* h2t; %生成下支路频带信号figure(2);subplot(2,4,1);plot(t,s_16qam1);axis(tx1,tx2,ty1,ty2);xlabel(时间(S);ylabel(幅度);title(上支路已调信号时域波形图);grid;f2,s_16qam1f=T2F(t,s_16qam1);figure(2);subplot(2,4,5);plot(f2,10*log10(abs(s_16qam1f).2/T);ax
31、is(fx1,fx2,fy1,fy2);xlabel(频率(Hz);ylabel(功率谱密度(dB/Hz);title(上支路已调信号功率谱图);grid;figure(2);subplot(2,4,2);plot(t,s_16qam2);axis(tx1,tx2,ty1,ty2);xlabel(时间(S);ylabel(幅度);title(下支路已调信号时域波形图);grid;f3,s_16qam2f=T2F(t,s_16qam2);figure(2);subplot(2,4,6);plot(f3,10*log10(abs(s_16qam2f).2/T);axis(fx1,fx2,fy1,f
32、y2);xlabel(频率(Hz);ylabel(功率谱密度(dB/Hz);title(下支路已调信号功率谱图);grid;s_16qam=s_16qam1+s_16qam2;figure(2);subplot(2,4,3);plot(t,s_16qam);axis(tx1,tx2,ty1,ty2);xlabel(时间(S);ylabel(幅度);title(已调信号时域波形图);grid;f4,s_16qamf=T2F(t,s_16qam);figure(2);subplot(2,4,7);plot(f4,10*log10(abs(s_16qamf).2/T);axis(fx1,fx2,fy
33、1,fy2);xlabel(频率(Hz);ylabel(功率谱密度(dB/Hz);title(已调信号信号功率谱图);grid; % 信道加入高斯白噪声进行接受解调% 产生高斯白噪声 m1=1;p1=-10;noise =wgn(m1,Lt,p1);% 接收信号y_16qam = s_16qam + noise;figure(2);subplot(2,4,4);plot(t,y_16QAM);axis(tx1,tx2,ty1,ty2);xlabel(时间(S);ylabel(幅度);title(接收信号时域波形图);grid;f5,y_16qamf=T2F(t,y_16qam);figure(
34、2);subplot(2,4,8);plot(f5,10*log10(abs(y_16qamf).2/T);axis(fx1,fx2,fy1,fy2);xlabel(频率(Hz);ylabel(功率谱密度(dB/Hz);title(接收信号功率谱图);grid;% 相干解调% 通过乘法器1r_16QAM11 = y_16qam .* h1t;%通过低通滤波器b,a=butter(3,0.1);r_16QAM11=filter(b,a,r_16QAM11);figure(3);subplot(2,2,2);plot(t,r_16QAM11)axis(tx1,tx2,ty1,ty2);xlabel
35、(时间(S);ylabel(幅度);title(上支路通过低通滤波器信号时域波形图);grid;%抽样判决dd111=r_16QAM11(2*fc*N_sample:4*fc*N_sample:end); %在每个码元中间抽样 dd1111=;for n=1:N/4 if dd111(n)1 dd1111(n)=3; elseif 0dd111(n) & dd111(n)1 dd1111(n)=1; elseif -1dd111(n) & dd111(n)1 dd2222(n)=3; elseif 0dd222(n) & dd222(n)1 dd2222(n)=1; elseif -1dd22
36、2(n) & dd222(n)0 dd2222(n)=-1; else dd2222(n)=-3; endenddd22222=sigexpand(dd2222,4*fc*N_sample); d_NRZ23=conv(dd22222,gt2);figure(4);subplot(2,2,3);plot(t,d_NRZ23(1:Lt)axis(tx1,tx2,ty1,ty2);xlabel(时间(S);ylabel(幅度);title(下支路抽样判决后信号时域波形图);grid;% 4-2电平转换dd111111=;for p=1:N/4 if dd1111(p)=3 dd111111(2*p-1)=1; dd111111(2*p)=1; elseif dd1111(p)=1 dd111111(2*p-1)=1; dd111111(2*p)=-1; elseif dd1111(p)=-1 dd111111(2*p-1)=-1; dd111111(2*p)=1; elseif dd1111(p)=-3 dd111111(2*p-1)=-1; dd111111(2*p