function [numoferr]=qbpskberr(A,fc,d,snr,N_sample,N,Ts,df)
% 求误码率
% -------------------系统仿真参数------------------------------
%clc;clear all;close all;A=1;fc=1;snr=10;N_sample=100;N=10;Ts=1;df=0.01;
% A: 载波振幅 fc: 载波频率(Hz)
% snr: 信噪比dB N_sample: 每个码元的的采样点数
% N: 二进制码元数 Ts: 码元宽度
% df: 频率分辨率
% -------------------输出(返回)参数----------------------------
% numoferr: 误码率
% panjue: 恢复的二进制代码1用1表示,0用-1表示
% desingal: 恢复的数字基带信号
% t: 时域采样时序
%----------------------生成调制信号
B=1/Ts; % 传码率
f_start=fc-B; % BPF的下截至频率
f_cutoff=fc+B; % BPF的上截至频率
fs=fc*N_sample; % 系统采样频率
ts=Ts/fs; % 系统采样间隔
t=0:ts:N/2*Ts-ts; % 时间轴
Lt=length(t); % 时间序列长度
a=0;
% 产生二进制信源
%d=round(rand(1,N)); %输入二进制代码
for i=1:N/2
d1(i)=d(i*2-1);
d2(i)=d(i*2);
end %串并转化
p0=zeros(1,N_sample);
p1=ones(1,N_sample);
PN=[];
PN1=[];
PN2=[];
for i=1:N
if d(i)==1;
PN=[PN p1];
else
PN=[PN p0];
end
end
% t1=0:ts/20:N/2*Ts-ts/20;
for i=1:N/2
if d1(i)==1
PN1=[PN1 p1];
else
PN1=[PN1 p0];
end
if d2(i)==1
PN2=[PN2 p1];
else
PN2=[PN2 p0];
end
end
d_sjx1=2*PN1-1; % 单/双后的NRZ(不归零)信号
d_sjx2=2*PN2-1;
%figure(1);
%subplot(2,2,1);
%plot(t,d_sjx1);
%axis([-0.1 5.3 -1.5 1.5]);
%ylabel('NRZ1');
%subplot(2,2,2);
%plot(t,d_sjx2);
%axis([-0.1 5.3 -1.5 1.5]);
%ylabel('NRZ2');
%----------------对数字基带信号进行4PSK调制-----------------
ht1=A*sin(2*pi*fc*t); % 载波
ht2=A*cos(2*pi*fc*t);
%subplot(2,2,3);
%plot(t,ht1);
%axis([-0.1 5.3 -1.5 1.5]);
%ylabel('载波1');
%subplot(2,2,4);
%plot(t,ht2);
%axis([-0.1 5.3 -1.5 1.5]);
%ylabel('载波2');
s_4psk=d_sjx1(1:Lt).*ht1+d_sjx2(1:Lt).*ht2; % 生成4PSK信号
%figure(2);
%subplot(2,2,1);
%plot(t,s_4psk);
%axis([-0.1 5.3 -1.5 1.5]);
%ylabel('s_4psk');
%----------------生成高斯白噪声噪声-------------------------
snr_lin=10^(snr/10); % 换算成倍数
signal_power=dot(s_4psk(1:Lt),s_4psk(1:Lt))/Lt; % 求出信号平均功率
noise_power=(signal_power*fc*N_sample)/(snr_lin*2*(f_cutoff-f_start)); %求出噪声方差
noise_std=sqrt(noise_power); % 求出噪声均方差
noise=noise_std.*randn(1,Lt); % 以噪声方均根作为幅度产生高斯白噪声
%subplot(2,2,2);
%plot(t,noise);
%axis([-0.1 5.3 -3.0 3.0]);
%ylabel('noise');
%---------------将已调信号送入信道-------------------------
r=s_4psk(1:Lt)+noise(1:Lt); % 叠加了噪声的已调信号
%subplot(2,2,3);
%plot(t,r);
%axis([-0.1 5.3 -5.0 5.0]);
%ylabel('已调信号');
[rf,r,df1,f]=T2F(r,ts,df,fs); % 接收信号的频谱
%subplot(2,2,4);
%plot(rf);
%axis([-10.1 10.3 -10.0 10.0]);
%ylabel('频谱序列');
%---------------在接收端先通过带通滤波器-------------------
[H,f]=bp_f(length(rf),f_start,f_cutoff,df1,fs,1); % 带通滤波器的频谱
DEM = H.*rf; % 带通滤波器输出信号的频谱
%figure(3);
%subplot(2,2,1);
%plot(DEM);
%axis([-10.1 10.3 -10.0 10.0]);
%ylabel('输出信号频谱序列');
[dem]=F2T(DEM,fs); % 滤波器的输出信号波形
%subplot(2,2,2)
%plot(dem);
%axis([-0.5 1024.5 -3.0 3.0]);
%ylabel('滤波器输出信号');
%---------------和本地载波相乘,进行相干解调---------------
der1=dem(1:Lt).*ht1; % 混频
der2=dem(1:Lt).*ht2; % 混频
%subplot(2,2,3);
%plot(der1);
%axis([-0.5 500.5 -3.0 3.0]);
%ylabel('解调输出信号1');
%subplot(2,2,4);
%plot(der2);
%axis([-0.5 500.5 -3.0 3.0]);
%ylabel('解调输出信号2');
[derf1,der1,df1,f]=T2F(der1,ts,df,fs); % 求混频后信号的功率谱
[derf2,der2,df2,f]=T2F(der2,ts,df,fs);
%figure(4);
%subplot(2,2,1);
%plot(derf1);
%axis([-10.5 10.5 -10.0 10.0]);
%ylabel('混频后信号1的功率谱');
%subplot(2,2,2);
%plot(derf2);
%axis([-10.5 10.5 -10.0 10.0]);
%ylabel('混频后信号2的功率谱');
%---------------再经过低通滤波器---------------------------
[LPF,f]=lp_f(length(derf1),B/2,df1,fs,1); % 产生LPF的传递函数
DM1 = LPF.*derf1; % 信号过理想低通滤波器输出的频谱
DM2 = LPF.*derf2; % 信号过理想低通滤波器输出的频谱
[dm1]=F2T(DM1,fs); % 理想低通滤波器的输出时域信号
[dm2]=F2T(DM2,fs); % 理想低通滤波器的输出时域信号
%subplot(2,2,3);
%plot(dm1);
%axis([-0.5 500.5 -1.3 1.3]);
%ylabel('过低通后的信号1');
%subplot(2,2,4);
%plot(dm2);
%axis([-0.5 500.5 -1.3 1.3]);
%ylabel('过低通后的信号2');
r_PN=[];
%---------------最后对LPF输出信号抽样判决------------------
panjue=zeros(1,N); % 建立存储判决值的矩阵
% 抽样判决,规则:大于0判1,小于0判-1
for i=1:N/2
n_t1=fc*N_sample*(i-1)+1;
n_t2=fc*N_sample*i;
n_t=fc*N_sample*(i-1)+fc*N_sample/2;
if (dm1(n_t))>=0; % 抽样判决时刻
panjue(2*i-1)=1; % 大于等于0判1
r_PN=[r_PN p1];
else
panjue(2*i-1)=-1; % 小于0判-1
r_PN=[r_PN p0];
end
if (dm2(n_t))>=0; % 抽样判决时刻
panjue(2*i)=1; % 大于等于0判1
r_PN=[r_PN p1];
else
panjue(2*i)=-1; % 小于0判-1
r_PN=[r_PN p0];
end
end
figure(5);
subplot(2,1,1);
plot(r_PN);
axis([0 1000 -1.1 1.1]);
ylabel('判决结果');
subplot(2,1,2);
plot(PN);
axis([0 1000 -1.1 1.1]);
ylabel('原始信号');
%----------------统计误码数-------------------------------
% plot(t,dm1(1:Lt),t,dm2(1:Lt),'.',t1,PN,t1,r_PN);
%numoferr=(sum(abs(panjue-d))/2)/N;
for i=1:N/2
if(panjue(2*i-1)==d1(i)&&panjue(2*i)==d2(i))
a=a+1;
end
end
numoferr=(N/2-a)/(N/2);