ofdm_hilbert.m

OFDM信道估计和均衡的仿真程序

%+----------------------------------------------------------------+
%|                                                                |
%|  Name: ofdm_Hilbert.m                                                 |
%|  Author: Mathew Liu, James Shen                                |
%|  Description: Basic OFDM model used for simulation purposes.   |    
%|  Using rcosine as transmit filter. Hilbert transform at        |
%|  recevier.                                                     |
%|                                                                |
%+----------------------------------------------------------------+
% We assume that with are using 64 length FFT with 64 subcarriers all carrying
% data. So there are no carriers specifically assigned to pilot.
clear all;
close all;
%=============== Some standard Hyperlan Params ==================
T = 50e-9;  % System baseband sampling period.  
fs = 1/T;   % System baseband sampling freq = 20MHz .
Tcp = 16*T; % CP period.
Tu = 64*T;  % Useful symbol period.
Ts = Tu+Tcp; % OFDM Symbol period 80 samples.
delta_f = 0.3125e6; % Frequency spacing.
time_scale = T:T:Ts;    % Time samples in baseband signals.
N = length(time_scale); % NUmber of time samples in baseband.
f1 = 1/Ts:1/Ts:N/(2*Ts);   % Frequencies in the 1st half of Nyquist circle.
f2 = -N/(2*Ts):1/Ts:-1/Ts; % Frequencies in the 2nd half of Nyquist circle.
freq_scale = [f1 f2];   % Baseband frequency scale.
fc = 2.4e9; % Carrier frequency
% If Rfs=2*fc then,
%       If fc=5.2, over=size(rt)/80 - 1=519;
%       If fc=2.4, over=239
% If Rfs=5*fc then,
%       If fc=5.2, over=size(rt)/80 - 1=1299;
%       If fc=2.4, over=599
Rfs = 5*fc; % Simulation frequency. Mimics real time behaviour for bandpass signal. 
            % Sampling at Nyquist rate of twice the highest freq.
RT = 1/Rfs;      % The real-time sampling period 
rt = RT:RT:Ts;   % The time samples in real-time bandpass
RN = length(rt); % Number of samples in real-time bandpass
f3 = 1/Ts:1/Ts:RN/(2*Ts);
f4 = -RN/(2*Ts):1/Ts:-1/Ts;
rf_scale = [f3 f4]; % Bandpass freq scale.
%============== Our params based on Hyperlan/2 =================
FFTLen = 64; % Length of FFT.
CPLen = 16;  % Length of Cyclic Prefix
M = 4;       % Bits encoded in a QAM symbol.

%+----------------------------------------------------------------+
%|                    OFDM   TRANSMITTER                          |
%+----------------------------------------------------------------+
%======== Data Generation ===============================
input = rand(1,FFTLen*M) > 0.5; % transmits one symbol

%======== Serial To Block (FFTLen, M) ============================
input_para = reshape(input,FFTLen,M);

%======== Baseband Modulation ===========================
input_symbols = qammod(input_para,FFTLen, M);
figure(1);
plot(input_symbols,'ob');grid;axis square;axis equal;

%================= IFFT =================================
input_time_para = ifft(input_symbols);

%================ Parallel to Serial ====================
input_time_serial = reshape(input_time_para,1,FFTLen); % Time domain signal

%=============  Gaurd Interval Insertion (GII) =================
input_ext = [input_time_serial((FFTLen-CPLen+1):FFTLen),input_time_serial];
figure(2);
subplot(2,1,1);
stem(time_scale, abs(input_ext)), title('Cycl. Ext. Signal (Time)');
subplot(2,1,2);
stem(freq_scale, abs(fft(input_ext))), title('Specturm of Cycl. Ext. Signal (Freq)');

%=================== Pulse shaping Tansmit Filter ===============================
%input_shaped = input_ext.*hamming(FFTLen+CPLen).';
over = 599; % Oversampling factor according to Simulation frequency.
input_over = kron(input_ext, [1, zeros(1,over)]);
Nover = length(input_over);
%pulse_shape = ones(1, over+1);
pulse_shape = rcosine(1,(over+1),'normal',0.35); % basic raised-cosine pulse-shape
[trash,pos] = max(pulse_shape); % raised cosine filter delay 'pos'
input_shaped_temp = conv(pulse_shape, input_over);

input_shaped = input_shaped_temp(pos:(pos+Nover-1));% Get rid of transient

%-----------------Over sampling scales---------------
overt = Ts/Nover:Ts/Nover:Ts;
f5 = 1/T:1/T:Nover/(2*T);
f6 = -Nover/(2*T):1/T:-1/T;
overf_scale = [f5 f6];
%================= Graphing Section ================================
figure(3);
subplot(2,1,1);
stem(overt(1:5000), abs(input_over(1:5000)), 'b'), title('Signal after over sampling'); hold on;
plot(overt(1:5000), abs(input_shaped(1:5000)), 'r'), title('Signal after pulse shaping.'); hold off;
subplot(2,1,2);
stem(overf_scale, abs(fft(input_shaped))), title('Spectrum of signal after pulse shaping');

%================== Upconverter to Carrier Freq ==========================

input_bandpass = 2*real(input_shaped.*exp(j*2*pi*fc*rt)); % Bandpass signal
figure(4);
subplot(2,1,1);
plot(rt(1:5000), input_bandpass(1:5000)), title('Bandpass signal (time)');
subplot(2,1,2);
stem(rf_scale, abs(fft(input_bandpass))), title('Spectrum of Bandpass Signal (Freq)');

%+----------------------------------------------------------------+
%|                    CHANNEL MODEL                               |
%+----------------------------------------------------------------+
%=============== Pass Through Rayleigh Channel ==========
output_ext = filter(1, 1, input_bandpass); % Does nothing at the moment.

%+----------------------------------------------------------------+
%|                    OFDM   RECEVIER                             |
%+----------------------------------------------------------------+
%================ Downconvert ==================================
rec_over_bad = output_ext;
figure(5);
subplot(2,1,1);
plot(rt(1:5000), abs(rec_over_bad(1:5000))), ...
    title('Received Baseband Signal (time)'), ...
    ylabel('Amplitude'), xlabel('Time (s)');
subplot(2,1,2);
rec_over_temp = abs(fft(rec_over_bad));
rec_over_temp = rec_over_temp./max(abs(rec_over_temp)); % Normalize spectrum for plotting.
stem(rf_scale, rec_over_temp), title('Spectrum of Baseband Signal (Freq)'), ...
    ylabel('Amplitude'), xlabel('Frequncy (Hz)'); hold on;

h = remez(200,[0.02 0.98],[1 1],'Hilbert'); % Design of a Hilbert filter of order 200
f2 = 0.5*([zeros(1,100) 1 zeros(1,100)]+ j*h); % Complex Phase splitter filter whose real part is a delay and
pos = 101; % imaginary part is a Hilbert transformer
F2 = freqz(f2,1,-pi:2*pi/RN:pi-2*pi/RN); % Frequency response of the complex filter
plot([f4, f3],abs(F2),'r'); hold off;

rec_over_temp = conv(f2,rec_over_bad); % Signal is filtered

rec_over = rec_over_temp(pos:(pos+Nover-1)); % Get rid of inital transient part.

rec_baseband = rec_over.*exp(-j*2*pi*fc*rt); % I/Q downconversion

figure(6);
subplot(2,1,1);
plot(rt(1:5000), abs(rec_baseband(1:5000))), title('Received Baseband Downconverted Signal (time)');
subplot(2,1,2);
stem(rf_scale, abs(fft(rec_baseband))), title('Spectrum of Baseband Downconverted Signal (Freq)'); hold on;

output_ext = rec_baseband(1:over+1:Nover);

%============== Gaurd Interval Removal (GIR) ==================
output_time_serial = output_ext((CPLen+1):FFTLen+CPLen);

%======== Serial To Parallel =================================
output_time_para = reshape(output_time_serial,FFTLen,1);

%================= FFT =================================
output_symbols = fft(output_time_para);

%================= Baseband Demodulation ===========================
figure(7);
plot(output_symbols,'ob');grid;axis square;axis equal;

output_para = qamdemod(output_symbols,FFTLen, M);

%================ Parallel to Serial ====================
output = reshape(output_para,1,FFTLen*M); % Time domain signal

%==================== BER Calculation =========================
errors = abs(input - output);
num_errors = sum(sum(errors));
BER = num_errors/(FFTLen*M)