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

主程序部分：

%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 ***************************

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

matlab自带了一个mimo信道的模块，在block simulink里，这个简帆旁短不了；

部分参考如下：

SNR=0%定义信噪比为0HZ

A=10^(SNR/10)%信噪比的单位转换关系式

Im=eye(1)%这里的"1"为发射天线和接收猛散天枝轿氏线两者数目少的根数

sum=0%初始化sum为0

for i=1:10000 %程序循环10000次

H=randn(1,4)%产生高斯信号源

Q=H*H'%求威沙特矩阵

sum=sum+log2(det(Im+A*Q/4))%10000次循环得到的信道容量总和

end

C=sum/10000 %10000次循环的平均值

---