message=randint(1,10000);%随机生成10000个0 1码
num=10;%重复实验次数
%...................................卷积编码...............................
G=[1 1 1 1 0 0 1;1 0 1 1 0 1 1];
k0=1;
L=size(G,2)/k0;
n0=size(G,1);
output=cnv_encd(G,k0,message); %卷积编码结果
% channel_output=output;
% %output1=output(1,1:length(output)-n0*(L-1));
% [decoder_output,survivor_state,cumulated_metric]=viterbi(G,k0,channel_output);
%.......将卷积编码序列转换为双极性序列.............................
for i=1:length(output)
if output(1,i)==0
output(1,i)=-1;
else
output(1,i)=1;
end
end
%.................MSK调制...................
data=output;
data_len=length(output); %码元个数
sample_number = 8; %采样个数
Rb = 1500000; %码元速率 采样时间20.6666s
fc = 6000000; %载波频率
%两个频率为6375000,5625000 带宽3750000
%................MSK基带调制...................%
[signal_out,I_out,Q_out] = mod_msk(data,data_len,sample_number,Rb);
%.............................................%
multi = fc/Rb;
I_temp=interp(I_out,multi);
Q_temp=interp(Q_out,multi);
Fs=Rb*sample_number*multi;
t=1/Fs:1/Fs:length(I_temp)*1/Fs;
signal_i=I_temp.*cos(2*pi*fc*t);
signal_q=Q_temp.*sin(2*pi*fc*t);
signal_mod=I_temp.*cos(2*pi*fc*t)-Q_temp.*sin(2*pi*fc*t);
N=sample_number*data_len*multi; %总采样点数
...............................................................%
SNR0=-10;
SNR1=10;
wumalvjuzhen1=zeros(1,SNR1-SNR0+1);
wumalvjuzhen1z=zeros(1,SNR1-SNR0+1);
for n=1:num; %.............................................................
for k=SNR0:SNR1
% snr=k;
% snr1=10^(snr/10); %将信噪比的值由dB转化为数值
signal_mod1=awgn(signal_mod,k,'measured');
N=300; % 滤波器的阶数为(N+1)
F=[0,fc-375000,fc+375000,Fs/2]*2/Fs;
A=[1,1,0,0];
lpf=firls(N,F,A);
[amp_lpf,w]=freqz(lpf);
I_dem=signal_mod1.*cos(2*pi*fc*t)*2;
I_dem=conv(I_dem,lpf);
I_dem=I_dem(N/2+1:N/2+length(I_temp));
Q_dem=signal_mod1.*sin(2*pi*fc*t)*2;
Q_dem=conv(Q_dem,lpf);
Q_dem=-Q_dem(N/2+1:N/2+length(Q_temp));
I_dem_out=zeros(1,length(I_dem)/multi); % 抽取
Q_dem_out=zeros(1,length(Q_dem)/multi);
for i=1:length(I_dem_out)
I_dem_out(i)=I_dem(multi*(i-1)+1);
Q_dem_out(i)=Q_dem(multi*(i-1)+1);
end;
%..................................................................
%差分解调
demod_data = zeros(1,data_len);
demod_data(1) = Q_dem_out(sample_number);
for i = 2:2:data_len
demod_data(i) = -I_dem_out(i*sample_number)*Q_dem_out((i-1)*sample_number);
end
for i = 3:2:data_len
demod_data(i) = Q_dem_out(i*sample_number)*I_dem_out((i-1)*sample_number);
end
demod_data=demod_data>0;
demod_data=2*demod_data-1;
%demod_data即为解调序列
for i=1:length(demod_data)
if demod_data(1,i)==-1
demod_data(1,i)=0;
else
demod_data(1,i)=1;
end
end
%将解调出来的信息序列进行维特比译码,decoer_output是译码输出,即应为message
[decoder_output,survivor_state,cumulated_metric]=viterbi(G,k0,demod_data);
% end
wumashu=length(find(decoder_output-message));
wumalv=wumashu/length(message); %计算误码率
wumalvjuzhen1(k-SNR0+1)=wumalv;
end
wumalvjuzhen1z=wumalvjuzhen1z+wumalvjuzhen1;
end;
wumalvjuzhen1z=wumalvjuzhen1z/num;
for i=1:length(message)
if message(1,i)==0
message(1,i)=-1;
else
message(1,i)=1;
end
end
%.................MSK调制...................
data=message;
data_len=length(message); %码元个数
sample_number = 8; %采样个数
Rb = 1500000; %码元速率 采样时间20.6666s
fc = 6000000; %载波频率
%两个频率为6375000,5625000 带宽3750000
%................MSK基带调制...................%
[signal_out,I_out,Q_out] = mod_msk(data,data_len,sample_number,Rb);
%.............................................%
multi = fc/Rb;
I_temp=interp(I_out,multi);
Q_temp=interp(Q_out,multi);
Fs=Rb*sample_number*multi;
t=1/Fs:1/Fs:length(I_temp)*1/Fs;
signal_i=I_temp.*cos(2*pi*fc*t);
signal_q=Q_temp.*sin(2*pi*fc*t);
signal_mod=I_temp.*cos(2*pi*fc*t)-Q_temp.*sin(2*pi*fc*t);
N=sample_number*data_len*multi; %总采样点数
...............................................................%
SNR0=-10;
SNR1=10;
wumalvjuzhen2=zeros(1,SNR1-SNR0+1);
wumalvjuzhen2z=zeros(1,SNR1-SNR0+1);
for i=1:length(message)
if message(1,i)==-1
message(1,i)=0;
else
message(1,i)=1;
end
end
for n=1:num
for k=SNR0:SNR1
% snr=k;
% snr1=10^(snr/10); %将信噪比的值由dB转化为数值
signal_mod1=awgn(signal_mod,k,'measured');
N=300; % 滤波器的阶数为(N+1)
F=[0,fc-375000,fc+375000,Fs/2]*2/Fs;
A=[1,1,0,0];
lpf=firls(N,F,A);
[amp_lpf,w]=freqz(lpf);
I_dem=signal_mod1.*cos(2*pi*fc*t)*2;
I_dem=conv(I_dem,lpf);
I_dem=I_dem(N/2+1:N/2+length(I_temp));
Q_dem=signal_mod1.*sin(2*pi*fc*t)*2;
Q_dem=conv(Q_dem,lpf);
Q_dem=-Q_dem(N/2+1:N/2+length(Q_temp));
I_dem_out=zeros(1,length(I_dem)/multi); % 抽取
Q_dem_out=zeros(1,length(Q_dem)/multi);
for i=1:length(I_dem_out)
I_dem_out(i)=I_dem(multi*(i-1)+1);
Q_dem_out(i)=Q_dem(multi*(i-1)+1);
end;
%..................................................................
%差分解调
demod_data = zeros(1,data_len);
demod_data(1) = Q_dem_out(sample_number);
for i = 2:2:data_len
demod_data(i) = -I_dem_out(i*sample_number)*Q_dem_out((i-1)*sample_number);
end
for i = 3:2:data_len
demod_data(i) = Q_dem_out(i*sample_number)*I_dem_out((i-1)*sample_number);
end
demod_data=demod_data>0;
demod_data=2*demod_data-1;
%demod_data即为解调序列
for i=1:length(demod_data)
if demod_data(1,i)==-1
demod_data(1,i)=0;
else
demod_data(1,i)=1;
end
end
% end
wumashu=length(find(demod_data-message));
wumalv=wumashu/length(message); %计算误码率
wumalvjuzhen2(k-SNR0+1)=wumalv;
end
wumalvjuzhen2z=wumalvjuzhen2z+wumalvjuzhen2;
end
wumalvjuzhen2z=wumalvjuzhen2z/num;
x=SNR0:SNR1;
figure;
semilogy(x,wumalvjuzhen1z,'r-*',x,wumalvjuzhen2z,'b-o');
axis([-10,8,0,1]);
%semilogy(x,wumalvjuzhen1,'k-*',x,wumalvjuzhen2,'b-*',x,wumalvjuzhen3,'g-*');
% % semilogy(x,wumalvjuzhen,'k',x,pe,'b');
% %semilogy(x,wumalvjuzhen1,'k-*',x,wumalvjuzhen2,'b-');
xlabel('SNR');
ylabel('误码率');
% legend('带宽最宽','一半','0.1');
title('卷积编码Viterbi译码误码率曲线');
legend('卷积编码','无信道编码');
%legend('1','0.1');
grid on;
% figure;
% subplot(3,1,1);hua_fft(ganrao1,Fs,1);
% subplot(3,1,2);hua_fft(ganrao2,Fs,1);
% subplot(3,1,3);hua_fft(ganrao3,Fs,1);
% figure;
% subplot(2,1,1);hua_fft(signal_mod,Fs,1);title('DS/MSK信号频谱');
% subplot(2,1,2);hua_fft(ganrao1,Fs,1);title('窄带干扰频谱');
%
% wumashu=length(find(demod_data-output));
% wumalv=wumashu/length(output); %计算误码率