ofdm_system.m.bak

自己做的一个OFDM系统的波形信道仿真,分成各种模块实现

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% for OFDMA, OFDMA baseband modulation, ifft, RF, channel distortion, RF demodulation, fft,
% decode, and spectral analysis

% author: 谢雨翔
% 提高0502班
% 学号：012005012722
%
% 注：
%		（1）这是主程序代码，点RUN即可
% 	（2）蒙特卡罗仿真需要一定的时间，在RUN的过程中可能需要耐心等待半分钟左右（视精度而定）
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all; close all;

% parameters:
Nc = 16;      %number of subcarriers
Ns = 6;    %number of consecutive symbols in time domain
df = 10;    %subcarrier spacing
% Ns_sample = 20;  %number of samples per symbol本来只需要 N_sample = N 的，但是高阶的看高频是多少倍低频段,为了显示模拟波形
fc = 400;
K = 1;      % K = 1: localized OFDM, K > 1: identically distributed OFDM
TransmitIter = 50;      % Transmitter transmit bits Interation...
fs = Nc * 80;
dt = 1 / fs;                %sampling time interval

A = 3;
EbN0 = 1 : 0.01 : 1.2;
sigma = 10 .^ (-EbN0);

Tu = 1 / df;         %time-domain: one symbol transmited, frequency-domain: N symbols transmitted modulated to N different subcarriers
t = 0 : dt : Ns * Tu - dt;
T = Ns * Tu;        % total time

% error_vector = [ ];
% errors = 0;
for n = 1 : length(sigma)
    Errors(n) = 0;
    
    for count = 1 : TransmitIter

        tx_d = (randn(1, Ns * Nc) > 0);
        %demultiplex information rate as 1/N initial rate:
        d = reshape(tx_d, Nc, Ns); %1 0 1 0 0 1 1 0 becomes as 1 0 1 0 _ 0 1 1 0 _ ... _ x x x x 
        % d = d';     %d becomes a Nc * Ns matrix, with each row corresponding to one out of four paths
        d = d * 2 - 1; %generate information digit as -1 or 1 % BPSK

        test_expand = Tu/dt /Nc;
        test_expand2 = Tu/dt;
        %滤波恢复 D/A 数字到模拟 滤波器
        sdt = Tu / Nc;
        t1 = -Tu : dt : Tu-dt; %length is 2 * Tu/dt
        gt = sinc(t1 / sdt);
        % stem(t1, gt);
        for k = 1 : Ns
            ofdm_symbol(:, k) = ifft(d(:, k), Nc);
            sI(:, k) = real(ofdm_symbol(:, k));     %时域，一个sI(:,k)表示一个Tu中的Nc个数据的时域形式，size 为 Nc，但是一个Tu中有Tu/dt个数据绘图样点
            sQ(:, k) = imag(ofdm_symbol(:, k));

            sI2(:, k) = sigexpand( sI(:, k), Tu/dt / Nc );
            sIt(:, k) = conv(sI2(:, k), gt );       % D/A output
            sQ2(:, k) = sigexpand( sQ(:, k), Tu/dt / Nc );
            sQt(:, k) = conv(sQ2(:, k), gt );       % D/A output

            sIt2(:, k) = sIt( length(t1)/2 + 1: length(t1)/2 + Tu/dt , k);%?
            sQt2(:, k) = sQt( length(t1)/2 + 1: length(t1)/2 + Tu/dt , k);%?
        end

        if count == 1 && n == 1          %in next times, no need to draw plot again to demonstration
            figure(1);
            subplot(411);
            stem(0 : dt : Tu-dt, sI2(1 : Tu/dt, 1) );grid on;
            title('同相分量数字域波形');
            subplot(412);
            plot(0 : dt : Tu-dt, sIt2(:, 1));grid on;
            title('同相分量模拟域波形');
            subplot(413);
            stem(0 : dt : Tu-dt, sQ2(1 : Tu/dt, 1) );grid on;
            title('正交分量数字域波形');
            subplot(414);
            plot(0 : dt : Tu-dt, sQt2(:, 1));grid on;
            title('正交分量模拟域波形');
        end

        sIt3 = reshape(sIt2, 1, Ns * test_expand2);
        sQt3 = reshape(sQt2, 1, Ns * test_expand2);

        if count == 1 && n == 1
            figure(2);
            subplot(211);
            plot(t, sIt3);
            title('发送Ns帧In-phase分量');
            subplot(212);
            plot(t, sQt3);
            title('发送Ns帧Qu-phase分量');
        end
        %高频调制
        st = sIt3 .* cos(2 * pi * fc * t) - sQt3 .* sin(2 * pi * fc * t);
        st = A * st  + sigma(n) * rand(1, length(t) ); %有噪声
%         st = A * st;                              %没有噪声（测试）

        if count == 1 && n == 1
            figure(3);
            plot(t, st);
            xlabel('t'); ylabel('RF波形');
        end

        %高频解调
        rIt = st .* cos(2 * pi * fc * t);
        rQt = -st .* sin(2 * pi * fc * t);
        rIt = rIt - mean(rIt);
        rQt = rQt - mean(rQt);

        [f, rIf] = T2F(t, rIt);
        [t, rIt] = lpf(f, rIf, 2 * Nc * df);
        [f, rQf] = T2F(t, rQt);
        [t, rQt] = lpf(f, rQf, 2 * Nc * df);

        if count == 1 && n == 1
            figure(4);
            subplot(211);
            plot(t, rIt);
            subplot(212);
            plot(t, rQt);
        end

        rt = rIt + j * rQt;
        for k = 1 : Ns
            rt2 = rt((k-1)*Tu/dt+1 : k*Tu/dt);
            t2 = t((k-1)*Tu/dt+1 : k*Tu/dt);
        %     r = rI + j*rQ;
            r = rt2(1 : Tu/dt/Nc : Tu/dt);
            rf = fft(r, Nc);
            dd1 = sign(real(rf));
        %     [f2, rf(:,k)] = T2F( t2, rt2(:,k) );
            dd(:,k) = dd1';
            dd(:,k) = (dd(:,k) + 1)/2;
        end

        if count == 1 && n == 1
            figure(5);
            subplot(211);
            plot(t2, real(rt2));
            axis([0 t2(end) -1 0.7]);
            subplot(212);
            plot(t2, imag(rt2));
            axis([0 t2(end) -1 0.7]);
        end

        rx_d = reshape(dd, 1, Ns * Nc);
        error_vector = mod(tx_d + rx_d, 2);
        Errors(n) = Errors(n) + sum(error_vector);
    end
    Errors(n) = Errors(n)/(TransmitIter * Nc * Ns);
end

figure(6);
semilogy(EbN0,Errors,'b*-');hold on;grid on;