Question

Matlab编写BPSK信号通过AWGN信道的仿真过程

求一份完整的Matlab编写BPSK信号通过AWGN信道仿真过程的程序。最好能有注释。在线求啊

Answer 1

主程序部分：
%programm 3-1

%bpsk.m
%
% Simulation program to realize BPSK transmission system
%
% Programmed by H.Harada and T.Yamamura,
%

%******************** Preparation part **********************

sr=256000.0; % Symbol rate
ml=1; % Number of modulation levels
br=sr.*ml; % Bit rate (=symbol rate in this case)
nd = 1000; % Number of symbols that simulates in each loop
ebn0=3; % Eb/N0
IPOINT=8; % Number of oversamples

%******************* Filter initialization ********************

irfn=21; % Number of filter 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

%******************** START CALCULATION *********************
nloop=100; % Number of simulation loops

noe = 0; % Number of error data
nod = 0; % Number of transmitted data

for iii=1:nloop

%******************** Data generation ***********************

data=rand(1,nd)>0.5; % rand: built in function

%******************** BPSK Modulation ***********************

data1=data.*2-1;
data2 = oversamp( data1, nd , IPOINT) ;
data3 = conv(data2,xh); % conv: built in function

%****************** Attenuation Calculation *****************

spow=sum(data3.*data3)/nd;
attn=0.5*spow*sr/br*10.^(-ebn0/10);
attn=sqrt(attn);

%********************** Fading channel **********************

% Generated data are fed into a fading simulator
% In the case of BPSK, only Ich data are fed into fading counter
% [ifade,qfade]=sefade(data3,zeros(1,length(data3)),itau,dlvl,th1,n0,itnd1,now1,length(data3),tstp,fd,flat);

% Updata fading counter
%itnd1 = itnd1+ itnd0;

%************ Add White Gaussian Noise (AWGN) ***************

inoise=randn(1,length(data3)).*attn; % randn: built in function
data4=data3+inoise;
data5=conv(data4,xh2); % conv: built in function

sample=irfn*IPOINT+1;
data6 = data5(sample:8:8*nd+sample-1);

%******************** BPSK Demodulation *********************

demodata=data6 > 0;

%******************** Bit Error Rate (BER) ******************

noe2=sum(abs(data-demodata)); % sum: built in function
nod2=length(data); % length: built in function
noe=noe+noe2;
nod=nod+nod2;

fprintf('%d\t%e\n',iii,noe2/nod2); %fprintf:built in function
end % for iii=1:nloop

%********************** Output result ***********************

ber = noe/nod;
fprintf('%d\t%d\t%d\t%e\n',ebn0,noe,nod,noe/nod);

%******************** end of file ***************************

下面是子程序：
% Program 3-2
% oversamp.m
%
% Insert zero data to input data
%
% Programmed by H.Harada
%

function [out] = oversamp( indata, nsymb , sample)

%****************** variables *************************
% indata : input sequence
% nsymb : Number of symbols
% sample : Number of oversample
% *****************************************************

out=zeros(1,nsymb*sample);
out(1:sample:1+sample*(nsymb-1))=indata;

%******************** end of file ***************************

子程序：
% Program 3-3
% hrollfcoef.m
%
% Generate coefficients of Nyquist filter
%
% programmed by H.Harada
%

function [xh] = hrollfcoef(irfn,ipoint,sr,alfs,ncc)

%****************** variables *************************
% irfn : Number of symbols to use filtering
% ipoint : Number of samples in one symbol
% sr : symbol rate
% alfs : rolloff coeficiense
% ncc : 1 -- transmitting filter 0 -- receiving filter
% *****************************************************

xi=zeros(1,irfn*ipoint+1);
xq=zeros(1,irfn*ipoint+1);

point = ipoint;
tr = sr ;
tstp = 1.0 ./ tr ./ ipoint;
n = ipoint .* irfn;
mid = ( n ./ 2 ) + 1;
sub1 = 4.0 .* alfs .* tr; % 4*alpha*R_s

for i = 1 : n

icon = i - mid;
ym = icon;

if icon == 0.0
xt = (1.0-alfs+4.0.*alfs./pi).* tr; % h(0)
else
sub2 =16.0.*alfs.*alfs.*ym.*ym./ipoint./ipoint;
if sub2 ~= 1.0
x1=sin(pi*(1.0-alfs)/ipoint*ym)./pi./(1.0-sub2)./ym./tstp;
x2=cos(pi*(1.0+alfs)/ipoint*ym)./pi.*sub1./(1.0-sub2);
xt = x1 + x2; % h(t) plot((1:length(xh)),xh)
else % (4alphaRst)^2 = 1plot((1:length(xh)),xh)
xt = alfs.*tr.*((1.0-2.0/pi).*cos(pi/4.0/alfs)+(1.0+2.0./pi).*sin(pi/4.0/alfs))./sqrt(2.0);
end % if sub2 ~= 1.0
end % if icon == 0.0

if ncc == 0 % in the case of receiver
xh( i ) = xt ./ ipoint ./ tr; % normalization
elseif ncc == 1 % in the case of transmitter
xh( i ) = xt ./ tr; % normalization
else
error('ncc error');
end % if ncc == 0

end % for i = 1 : n

%******************** end of file ***************************

就是这三个了，不知道能不能帮助你、、、