clear all;
clc;
Ns=8;%做3次小波包分解,得到8路信号
para=128;%小波包每一路数据符号的长度
N=1024;%小波包逆变换后数据的长度
bit=2;%每个符号的bit数
signal=round(rand(1,para*Ns*bit));%源信号
%调制和串并变换
din2=1-2*signal;
din_temp=reshape(din2,2,length(signal)/2);
for b=1:length(signal)/2
QPSKdata(b)=din_temp(1,b)+j*din_temp(2,b);
end
complex_qam_data=reshape(QPSKdata,8,length(QPSKdata)/8);
x=complex_qam_data;
%////////////////////调制和串并变换的实现////////////////////////////
%////////////////////显示调制后的星座图///////////////////////////////
figure(1);
plot(complex_qam_data,'*r');%16qam调制后星座图
axis([-4,4,-4,4]);
grid on
%////////////////////显示调制后的星座图///////////////////////////////
%//////////////////////wprec模块//////////////////////////////////////
%/////////////////////第三级重构/////////////////////////////////////
[ld,hd]=wfilters('db1','d');%求得db4小波的低通和高通分析滤波器
t30=x(1,:);
t31=x(2,:);
t32=x(3,:);
t33=x(4,:);
t34=x(5,:);
t35=x(6,:);
t36=x(7,:);
t37=x(8,:);
ld3=[ld,zeros(1,N/4-length(ld))];%第三级分解低通滤波器
hd3=[hd,zeros(1,N/4-length(hd))];%第三级分解高通滤波器
lr3=ld3(end:-1:1);
hr3=hd3(end:-1:1);
lr3=circshift(lr3',1)';
hr3=circshift(hr3',1)';
t30=dyadup(t30);
t31=dyadup(t31);
t32=dyadup(t32);
t33=dyadup(t33);
t34=dyadup(t34);
t35=dyadup(t35);
t36=dyadup(t36);
t37=dyadup(t37);
t30=t30(1:length(t30)-1);
t31=t31(1:length(t31)-1);
t32=t32(1:length(t32)-1);
t33=t33(1:length(t33)-1);
t34=t34(1:length(t34)-1);
t35=t35(1:length(t35)-1);
t36=t36(1:length(t36)-1);
t37=t37(1:length(t37)-1);
t30=ifft(fft(lr3).*fft(t30));
t31=ifft(fft(hr3).*fft(t31));
t20=t30+t31;
t32=ifft(fft(lr3).*fft(t32));
t33=ifft(fft(hr3).*fft(t33));
t21=t32+t33;
t34=ifft(fft(lr3).*fft(t34));
t35=ifft(fft(hr3).*fft(t35));
t22=t34+t35;
t36=ifft(fft(lr3).*fft(t36));
t37=ifft(fft(hr3).*fft(t37));
t23=t36+t37;
%////////////////////第二级重构/////////////////////////////////////
ld2=[ld,zeros(1,N/2-length(ld))];
hd2=[hd,zeros(1,N/2-length(hd))];
t20=dyadup(t20);
t21=dyadup(t21);
t22=dyadup(t22);
t23=dyadup(t23);
t20=t20(1:length(t20)-1);
t21=t21(1:length(t21)-1);
t22=t22(1:length(t22)-1);
t23=t23(1:length(t23)-1);
lr2=ld2(end:-1:1);
hr2=hd2(end:-1:1);
lr2=circshift(lr2',1)';
hr2=circshift(hr2',1)';
t20=ifft(fft(lr2).*fft(t20));
t21=ifft(fft(hr2).*fft(t21));
t10=t20+t21;
t22=ifft(fft(lr2).*fft(t22));
t23=ifft(fft(hr2).*fft(t23));
t11=t22+t23;
%///////////////第一级重构///////////////////////////////////////////
ld1=[ld,zeros(1,N-length(ld))];
hd1=[hd,zeros(1,N-length(hd))];
t10=dyadup(t10);
t11=dyadup(t11);
t10=t10(1:length(t10)-1);
t11=t11(1:length(t11)-1);
lr1=ld1(end:-1:1);
hr1=hd1(end:-1:1);
lr1=circshift(lr1',1)';
hr1=circshift(hr1',1)';
t10=ifft(fft(lr1).*fft(t10));
t11=ifft(fft(hr1).*fft(t11));
sig=t10+t11;%在信道中传输的信号为sig
%/////////////////////wprec模块///////////////////////////////////////
%%/////////////////添加噪声///////////////////////////////////////////
Tx_signal_power=var(sig);%发送信号功率
a=[0:0.2:50];
for b=1:length(a)
SNR=a(b);
linear_SNR=10^(SNR/10);%线性信噪比
noise_sigma=Tx_signal_power/linear_SNR;
noise_scale_factor=sqrt(noise_sigma);%标准差sigma
noise=randn(1,1024)*noise_scale_factor;%产生正态分布噪声序列
%noise=wgn(1,length(windowed_Tx_data),noise_sigma,'complex');%产生复高斯白噪声信号
Rx_data=sig+noise;%接收到的信号加噪声
y=Rx_data;
%/////////////////添加噪声///////////////////////////////////////////
%////////////////wpdec模块////////////////////////////////////////////
%%////////////第一级分解/////////////////////////
t10=ifft(fft(y).*fft(ld1));
t11=ifft(fft(y).*fft(hd1));
t10=dyaddown(t10);
t11=dyaddown(t11);
%%%%%%%%%%%第一级分解////////////////////////////
%//////////////////第二级分解///////////////////////
t20=ifft(fft(t10).*fft(ld2));
t21=ifft(fft(t10).*fft(hd2));
t22=ifft(fft(t11).*fft(ld2));
t23=ifft(fft(t11).*fft(hd2));
t20=dyaddown(t20);
t21=dyaddown(t21);
t22=dyaddown(t22);
t23=dyaddown(t23);
%///////////////第二级分解////////////////////////////
%//////////////第三级分解////////////////////////////
t30=ifft(fft(t20).*fft(ld3));
t31=ifft(fft(t20).*fft(hd3));
t32=ifft(fft(t21).*fft(ld3));
t33=ifft(fft(t21).*fft(hd3));
t34=ifft(fft(t22).*fft(ld3));
t35=ifft(fft(t22).*fft(hd3));
t36=ifft(fft(t23).*fft(ld3));
t37=ifft(fft(t23).*fft(hd3));
t30=dyaddown(t30);
t31=dyaddown(t31);
t32=dyaddown(t32);
t33=dyaddown(t33);
t34=dyaddown(t34);
t35=dyaddown(t35);
t36=dyaddown(t36);
t37=dyaddown(t37);
Y1=[t30;t31;t32;t33;t34;t35;t36;t37];
Rx_phase=angle(Y1);
Rx_mag=abs(Y1)
%figure;
%polar(Rx_phase,Rx_mag,'bd');%极坐标坐标下面画出接收信号的星座图
%////////////////wpdec模块////////////////////////////////////////////
%//////////////////接收信号的星座图///////////////////////////////////////
[M,N]=pol2cart(Rx_phase,Rx_mag);
Rx_complex_carrier_matrix=complex(M,N);
%figure;
%plot(Rx_complex_carrier_matrix,'*r');%XY坐标接收信号的星座图
%axis([-4 4 -4 4]);
%grid on;
%//////////////////接收信号的星座图///////////////////////////////////////
%//////////////////////QPSK解调//////////////////////////
Y=reshape(Y1,1,1024);
y(1,:)=real(Y);
y(2,:)=imag(Y);
yout=reshape(y,1,2048);
yout=-1.*yout;
ReSig=yout>0.001;
%//////////////////////QPSK解调//////////////////////////
%//////////////////计算误码率///////////////////////////////////////////
Signal_Rx=ReSig;
%figure(9);
%subplot(211);
%stem(signal(1:100));
% subplot(212);
% stem(Signal_Rx(1:100));
bit_errors=find(Signal_Rx-signal);
bit_error_count=size(bit_errors,2);
ber(b)=bit_error_count/(para*Ns*bit);
end
figure;
plot(ber);
title('基于小波包变换的OFDM系统的误码率');
