ber.m

ofdm系统的信道估计算法

function p=BER(SNRindB)
%clear all;
%close all;
%fprintf('OFDM仿真\n\n');
IFFT_bin_length=1024;      %FFT变换长度
carrier_count=200;        %子载波数
bits_per_symbol=2;        %每符号的位数
symbols_per_carrier=50;   %每个载波所含的符号数
%SNR=input('SNR=');
SNR=SNRindB;
baseband_out_length=carrier_count*symbols_per_carrier*bits_per_symbol;   %200*50*2=20000
carriers=(1:carrier_count)+(floor(IFFT_bin_length/4)-floor(carrier_count/2)); %(1:200)+(256-100) ??????（157--356）
conjugate_carriers=IFFT_bin_length-carriers+2;     %??????（770--869）
%信号发射
baseband_out=round(rand(1,baseband_out_length));  %(1,20000)
convert_matrix=reshape(baseband_out,bits_per_symbol,length(baseband_out)/bits_per_symbol);  %(2,10000)
for k=1:(length(baseband_out)/bits_per_symbol)  %%(1:10000)
    modulo_baseband(k)=0;
    for i=1:bits_per_symbol  %%(1:2)
        modulo_baseband(k)=modulo_baseband(k)+convert_matrix(i,k)*2^(bits_per_symbol-i);   %%第一行的2倍加第二行，2(1)+(2)
    end
end
carrier_matrix=reshape(modulo_baseband,carrier_count,symbols_per_carrier)';  %%(200*50)'=(50*200)
%QDPSK调制
carrier_matrix=[zeros(1,carrier_count);carrier_matrix];  %%(51*200)
for i=2:(symbols_per_carrier+1)  % 2:51
    carrier_matrix(i,:)=rem(carrier_matrix(i,:)+carrier_matrix(i-1,:),2^bits_per_symbol);  %%由0，1，2，3组成的
end
carrier_matrix=carrier_matrix*((2*pi)/(2^bits_per_symbol)); %%乘以pi/2
[X,Y]=pol2cart(carrier_matrix,ones(size(carrier_matrix,1),size(carrier_matrix,2)));  %%扩展为51*400，前半部分X是后半部分Y的反。
complex_carrier_matrix=complex(X,Y);  %%综合51*200:X+i*Y
%加训练序列
training_symbols=[1,j,j,1,-1,1,j,j,1,-j,j,-1,-1,-1,j,-j,1,j,j,1,-1,1,j,j,1,-j,j,-1,-1,-1,j,-j,1,j,j,1,-1,1,j,j,1,-j,j,-1,-1,-1,j,-j,1,j,j,1,-1,1,j,j,1,-j,j,-1,-1,-1,j,-j,1,j,j,1,-1,1,j,j,1,-j,j,-1,-1,-1,j,-j,1,j,j,1,-1,1,j,j,1,-j,j,-1,-1,-1,j,-j,1,j,j,1,-1,1,j,j,1,-j,j,-1,-1,-1,j,-j,1,j,j,1,-1,1,j,j,1,-j,j,-1,-1,-1,j,-j,1,j,j,1,-1,1,j,j,1,-j,j,-1,-1,-1,j,-j,1,j,j,1,-1,1,j,j,1,-j,j,-1,-1,-1,j,-j,1,j,j,1,-1,1,j,j,1,-j,j,-1,-1,-1,j,-j,1,j,j,1,-1,1,j,j,1,-j,j,-1,-1,-1,j,-j,1,j,j,1,-1,1,j,j];                     %%200个
training_symbols=cat(1,training_symbols,training_symbols);  %%将训练序列重复写，扩展为两行
training_symbols=cat(1,training_symbols,training_symbols);  %%重写，扩展为四行（4*200）
complex_carrier_matrix=cat(1, training_symbols, complex_carrier_matrix); %%将训练序列加在调制信号上面(55*200)
IFFT_modulation=zeros(4+symbols_per_carrier+1,IFFT_bin_length);   %%(55*1024)=0
IFFT_modulation(:,carriers) = complex_carrier_matrix;   %%前200列赋值为加训练序列的调制信号
IFFT_modulation(:,conjugate_carriers)=conj(complex_carrier_matrix);   %%770--869列为其共轭
time_wave_matrix=ifft(IFFT_modulation');  %%%转置并进行IFFT变换
time_wave_matrix = time_wave_matrix';  %%再进行转置
for i = 1:4+symbols_per_carrier+1      %%%1：55
    windowed_time_wave_matrix(i,:)=real(time_wave_matrix(i,:));  %%取实部55*1024
end
ofdm_modulation=reshape(windowed_time_wave_matrix',1,IFFT_bin_length*(4+symbols_per_carrier+1));  %%%并转串1*（1024*55）
Tx_data=ofdm_modulation;
%信道
d1=4;a1=0.2;d2=5;a2=0.3;d3=6;a3=0.4;d4=7;a4=0.5;
copy1=zeros(size(Tx_data)); %%%1*（1024*55）
for i=1+d1:length(Tx_data)   %%5：（1024*55）
    copy1(i)=a1*Tx_data(i-d1);%%1：（1024*55-5）Tx_data后移4列*0.2
end
copy2=zeros(size(Tx_data));
for i=1+d2:length(Tx_data)
    copy1(i)=a2*Tx_data(i-d2);
end
copy3=zeros(size(Tx_data));
for i=1+d3:length(Tx_data)
    copy1(i)=a3*Tx_data(i-d3);
end
copy4=zeros(size(Tx_data));
for i=1+d4:length(Tx_data)
    copy1(i)=a4*Tx_data(i-d4);
end
Tx_data=Tx_data+copy1+copy2;   %%%延迟？

Tx_signal_power=var(Tx_data);  %%方差
linear_SNR=10^(SNR/10);   %%线性信噪比
noise_sigma =Tx_signal_power/linear_SNR;   %%%sigma=p^2/(S/N)
noise_scale_factor=sqrt(noise_sigma);      %%%scale_factor=sqrt(p^2/(S/N))
noise=randn(1,length(Tx_data))*noise_scale_factor;%%%高斯随机过程
Rx_Data=Tx_data+noise;   %%叠加噪声的接收信号
%信号接收
Rx_Data_matrix=reshape(Rx_Data,IFFT_bin_length,4+symbols_per_carrier+1);  %%%1024*55
Rx_spectrum=fft(Rx_Data_matrix);  %%FFT变换
Rx_carriers=Rx_spectrum(carriers,:)';   %%（157--356）
Rx_training_symbols=Rx_carriers((1: 4),:);   %%((1: 4),(157--356))
Rx_carriers=Rx_carriers((5:55),:);
% 信道估计
Rx_training_symbols=Rx_training_symbols./training_symbols;   %%%(4*200)/(1*200)
Rx_training_symbols_deno=Rx_training_symbols.^2;
Rx_training_symbols_deno=Rx_training_symbols_deno (1,:)+Rx_training_symbols_deno(2,:)+Rx_training_symbols_deno(3,:)+Rx_training_symbols_deno(4,:);
Rx_training_symbols_nume=Rx_training_symbols(1,:)+Rx_training_symbols(2,:)+Rx_training_symbols(3,:)+Rx_training_symbols(4,:);
Rx_training_symbols_nume =conj (Rx_training_symbols_nume);
Rx_training_symbols = Rx_training_symbols_nume/Rx_training_symbols_deno;
Rx_training_symbols_2 = cat(1, Rx_training_symbols,Rx_training_symbols);
Rx_training_symbols_4 = cat(1, Rx_training_symbols_2,Rx_training_symbols_2);
Rx_training_symbols_8 = cat(1, Rx_training_symbols_4,Rx_training_symbols_4);
Rx_training_symbols_16= cat(1, Rx_training_symbols_8, Rx_training_symbols_8);
Rx_training_symbols_32= cat(1, Rx_training_symbols_16, Rx_training_symbols_16);
Rx_training_symbols_48= cat(1, Rx_training_symbols_32, Rx_training_symbols_16);
Rx_training_symbols_50= cat(1, Rx_training_symbols_48, Rx_training_symbols_2) ;
Rx_training_symbols = cat(1, Rx_training_symbols_50,Rx_training_symbols);
%Rx_carriers=Rx_training_symbols.*Rx_carriers;
Rx_phase = angle (Rx_carriers)*(180 /pi);
phase_negative = find(Rx_phase<0);  
Rx_phase(phase_negative)=rem(Rx_phase(phase_negative)+360,360);
Rx_decoded_phase= diff(Rx_phase);
phase_negative=find(Rx_decoded_phase<0);
Rx_decoded_phase(phase_negative)=rem(Rx_decoded_phase(phase_negative)+360,360);
% QDPSK解调
base_phase=360/2^bits_per_symbol;
delta_phase = base_phase/2;
Rx_decoded_symbols = zeros(size(Rx_decoded_phase,1),size(Rx_decoded_phase,2));
%信号输出
for i= 1:(2^bits_per_symbol-1)
    center_phase = base_phase*i;
    plus_deha = center_phase + delta_phase;
    minus_deha = center_phase - delta_phase;
    decoded = find((Rx_decoded_phase <=plus_deha)&(Rx_decoded_phase>minus_deha));
    Rx_decoded_symbols(decoded)=i;
end
Rx_serial_symbols = reshape(Rx_decoded_symbols',1,size(Rx_decoded_symbols,1)*size(Rx_decoded_symbols,2));
for i= bits_per_symbol:-1:1
    if i~=1
        Rx_binary_matrix(i,:)=rem(Rx_serial_symbols,2);
        Rx_serial_symbols=floor(Rx_serial_symbols/2);
    else
        Rx_binary_matrix(i,:)=Rx_serial_symbols;
    end
end
baseband_in=reshape(Rx_binary_matrix,1,size(Rx_binary_matrix,1)*size(Rx_binary_matrix,2));
%误码率计算
bit_errors=find(baseband_in~=baseband_out);
bit_error_count=size(bit_errors,2);
total_bits=size(baseband_out,2);
%bit_error_rate = bit_error_count/total_bits;
p = bit_error_count/total_bits;
%fprintf('%f\n',bit_error_rate);
%图