matlab_带瑞利衰落信道_高斯白噪声的QPSK调制解调信道，可以得到功率谱波形和星座图

%******************** Preparation part ************************************* clear all; echo off; sr=256000.0; % Symbol rate ml=2; % ml:Number of modulation levels (BPSK:ml=1, QPSK:ml=2, 16QAM:ml=4) br=sr .* ml; % Bit rate nd = 100; % Number of symbols that simulates in each loop T=nd/sr; IPOINT=8; % Number of oversamples %************************* Filter initialization *************************** irfn=21; % Number of taps alfs=0.5; % Rolloff factor [xh] = hrollfcoef(irfn,IPOINT,sr,alfs,1); %Transmitter filter coefficients [xh2] = hrollfcoef(irfn,IPOINT,sr,alfs,0); %Receiver filter coefficients %******************* Fading initialization ******************** tstp=1/sr/IPOINT; % Time resolution itau = [0]; % Arrival time for each multipath normalized by tstp % direct wave. dlvl = [0]; % Mean power for each multipath normalized by direct wave.% direct wave. n0=[6]; % Number of waves to generate fading for each multipath. th1=[0.0]; % Initial Phase of delayed wave itnd0=nd*IPOINT*100; % Number of fading counter to skip itnd1=[1000]; % Initial value of fading counter now1=1; % Number of directwave + Number of delayed wave fd=160; % Maximum Doppler frequency [Hz] flat =1; % flat : flat fading or not %******************** START CALCULATION ************************************* ebn0=0:1:10;% Eb/N0 for tt=1:length(ebn0) nloop=500; % Number of simulation loops noe = 0; % Number of error data nod = 0; % Number of transmitted data for iii=1:nloop %*************************** Data generation ******************************** data1=rand(1,nd*ml)>0.5; % rand: built in function %*************************** QPSK Modulation ******************************** [ich,qch]=qpskmod(data1,1,nd,ml); [ich1,qch1]= compoversamp(ich,qch,length(ich),IPOINT); [ich2,qch2]= compconv(ich1,qch1,xh); ich20 = ich2(irfn*IPOINT/2+1:nd*IPOINT+irfn*IPOINT/2); ICH2 = fftshift(fft(ich20)); P2 = ICH2.*conj(ICH2)/T; %功率计算 %**************************** Attenuation Calculation *********************** spow=sum(ich2.*ich2+qch2.*qch2)/nd; % sum: built in function attn=0.5*spow*sr/br*10.^(-ebn0/10); attn=sqrt(attn); % sqrt: built in function %********************** Fading channel ********************** % Generated data are fed into a fading simulator [ifade,qfade]=sefade(ich2,qch2,itau,dlvl,th1,n0,itnd1,now1,length(ich2),tstp,fd,flat); % Updata fading counter itnd1 = itnd1+ itnd0; %********************* Add White Gaussian Noise (AWGN) ********************** [ich3,qch3]= comb(ich2,qch2,attn(tt)); % 当需要计算衰落时，ich2,qch2改为ifade,qfade [ich4,qch4]= compconv(ich3,qch3,xh2); syncpoint=irfn*IPOINT+1; ich5=ich4(syncpoint:IPOINT:length(ich4)); qch5=qch4(syncpoint:IPOINT:length(qch4)); ich40 = ich4(irfn*IPOINT/2+1:nd*IPOINT+irfn*IPOINT/2); ICH4 = fftshift(fft(ich40)); P4 = ICH4.*conj(ICH4)/T; %功率计算 %**************************** QPSK Demodulation ***************************** [demodata]=qpskdemod(ich5,qch5,1,nd,ml); %************************** Bit Error Rate (BER) **************************** noe2=sum(abs(data1-demodata)); % sum: built in function nod2=length(data1); % length: built in function noe=noe+noe2; nod=nod+nod2; end % for iii=1:nloop %********************** Output result *************************** ber(tt) = noe/nod; fprintf('%e\n',noe/nod); % fprintf: built in function end %******************** 计算理论值 *************************** % SNRindB2=0:0.1:30; % for i=1:length(SNRindB2), % snr1=exp(SNRindB2(i)*log(10)/10); % snr2=sqrt(1+1/snr1); % snr3=1-1/snr2; % theo(i)=0.5*snr3; % end; semilogy(ebn0,ber,'g*'); % hold on; % semilogy(SNRindB2,theo); hold on; SNRindB3=0:0.2:10; for i=1:length(SNRindB3), SNR=exp(SNRindB3(i)*log(10)/10); % signal to noise ratio snr=sqrt(SNR) theo_err_prb(i)=0.5*erfc(snr); % theoretical bit error rat end; semilogy(SNRindB3,theo_err_prb); % Hs=spectrum.welch; % figure(2) % subplot(2,1,1) %基带波形 % stem(data1); % subplot(2,1,2) % psd(Hs,data1,'Fs',tstp); %基带功率谱 % % figure(3) % subplot(2,1,1) %调制波形 % plot(ich2); % subplot(2,1,2) % psd(Hs,ich2,'Fs',tstp); %调制后功率谱 % % % figure(4) % subplot(2,1,1) % plot(ich4); % subplot(2,1,2) % psd(Hs,ich4,'Fs',tstp); %加噪后功率谱 % % scatterplot([(ich)',(qch)']); % % scatterplot([(ich5)',(qch5)']); %星座图