itu_channel.m

OFDMA 物理层开发的matlab 源码.飞思卡尔提供.对物理层开发的工程师有帮助!

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Property of Freescale
%  Freescale Confidential Proprietary
%  Freescale Copyright (C) 2005 All rights reserved
%  ------------------------------------------------------------------------
%  $RCSfile: itu_channel.m.rca $
%  $Date: Sun Oct 22 03:19:51 2006 $
%  Target: Matlab
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 802.16-2004 OFDMA PHY - ITU-R M.1225 Wideband channel model
%
%    Description: y=itu_channel(x,fc,fs,model,v,samp_index) implements the
%                 terrestial mobile RF channel model defined in ITU-R M.1225.  
% 
%    Inputs:
%           x => Channel input matrix. Each column is a continuous stream
%                of information.
%           fc => Carrier frequency in Hz.
%           fs => Sampling frequency in Hz.
%           model => 'officeA', 'officeB', 'officeFlat', 'pedestrianA', 
%                    'pedestrianB','pedestrianFlat','vehicularA', 
%                    'vehicularB', 'vehicularFlat'
%           v => Mobile velocity in km/h. Optional parameter; defaults to 
%                0, 3, & 120 km/h for office, pedestrian, & vehicular models.
%           samp_index => Vector containing sample time index for the first 
%                         sample of each column. This input defaults to zero
%                         if its not present.
%
%    Outputs:    y => Each column of y is the output of the channel for the
%                corresponding column in x.
%                 
%    Data Files:
%                N/A
%                   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [y, impulse_res] = itu_channel(x,fc,fs,model,v,samp_index)

if nargin<6
    samp_index=0;
end

switch model
    case 'officeA'
        d=[0 50 110 170 290 310]*1e-9;
        g=10.^([0 -3 -10 -18 -26 -32]/20);
        g=g/sqrt(sum(g.^2));
        if nargin==4
            v=0;
        end
    case 'officeB'
        d=[0 100 200 300 500 700]*1e-9;
        g=10.^([0 -3.6 -7.2 -10.8 -18.0 -25.2]/20);
        g=g/sqrt(sum(g.^2));
        if nargin==4
            v=0;
        end
    case 'officeFlat'
        d=[0]*1e-9;
        g=10.^([0]/20);
        g=g/sqrt(sum(g.^2));
        if nargin==4
            v=0;
        end
    case 'pedestrianA'
        d=[0 110 190 410]*1e-9;
        g=10.^([0 -9.7 -19.2 -22.8]/20);
        g=g/sqrt(sum(g.^2));
        if nargin==4
            v=3;
        end
    case 'pedestrianB'
        d=[0 200 800 1200 2300 3700]*1e-9;
        g=10.^([0 -0.9 -4.9 -8.0 -7.8 -23.9]/20);
        g=g/sqrt(sum(g.^2));
        if nargin==4
            v=3;
        end
    case 'pedestrianC'
        d =[0 666.66 1333.33 2000 2666.66 3333.33 4000 4666.66 5333.33]*1e-9;
        g =10.^([0 -2 -4 -8 -15 -16 -23 -24 -25]/20);
        g=g/sqrt(sum(g.^2));
    case 'pedestrianFlat'
        d=[0]*1e-9;
        g=10.^([0]/20);
        g=g/sqrt(sum(g.^2));
        if nargin==4
            v=3;
        end
    case 'vehicularA'
        d=[0 310 710 1090 1730 2510]*1e-9;
        g=10.^([0 -1.0 -9.0 -10.0 -15.0 -20.0]/20);
        g=g/sqrt(sum(g.^2));
        if nargin==4
            v=120;
        end
    case 'vehicularB'
        d=[0 300 8900 12900 17100 20000]*1e-9;
        g=10.^([-2.5 0 -12.8 -10.0 -25.2 -16.0]/20);
        g=g/sqrt(sum(g.^2));
        if nargin==4
            v=120;
        end
    case 'vehicularFlat'
        d=[0]*1e-9;
        g=10.^([0]/20);
        g=g/sqrt(sum(g.^2));
        if nargin==4
            v=120;
        end
     case 'SUI-1'
        d=[0 0.4 0.8]*1e-6;
        g=10.^([0 -15 -20]/20);
        g=g/sqrt(sum(g.^2));
        if nargin==4
            v=120;
        end
     case 'SUI-2'
        d=[0 0.5 1]*1e-6;
        g=10.^([0 -12 -15]/20);
        g=g/sqrt(sum(g.^2));
        if nargin==4
            v=120;
        end
     case 'SUI-3'
        d=[0 0.5 1]*1e-6;
        g=10.^([0 -5 -10]/20);
        g=g/sqrt(sum(g.^2));
        if nargin==4
            v=120;
        end
     case 'SUI-4'
        d=[0 2 4]*1e-6;
        g=10.^([0 -4 -8]/20);
        g=g/sqrt(sum(g.^2));
        if nargin==4
            v=120;
        end
     case 'SUI-5'
        d=[0 5 10]*1e-6;
        g=10.^([0 -5 -10]/20);
        g=g/sqrt(sum(g.^2));
        if nargin==4
            v=120;
        end
     case 'SUI-6'
        d=[0 14 20]*1e-6;
        g=10.^([0 -10 -14]/20);
        g=g/sqrt(sum(g.^2));
        if nargin==4
            v=120;
        end
        
    otherwise
        disp('Unknown Channel Model')
        y=[];
        return
end

% Convert speed from km/h to m/s
v=v/3.6;

% Calculate doppler shift
fd=v*fc/3e8;

% Double the sample rate of the input signal. This prevents the
% delay.m function used in rayleigh.m from introducing distortion
% by attenuating the highest frequencies in x.

[InRows InColumns]=size(x);
delta = [1; zeros(InRows-1,1)]; 

total_rayleigh = rayleigh(1/fs,fd,InRows*InColumns,length(d));

% special case for k=1
% first symbol is without ISI
current_rayleigh= total_rayleigh((1:InRows),:);
y(:,1) = sum((current_rayleigh.*delay(x(:,1),d,fs))*g',2);
impulse_res(:,1) = sum((current_rayleigh.*delay(delta,d,fs))*g',2);

for k=2:InColumns,
    current_rayleigh1 = total_rayleigh((InRows*(k-1)+1:InRows*k),:);
    current_rayleigh2 = total_rayleigh((InRows*(k-2)+1:InRows*k),:);
    prevSym = x(:,k-1);
    currSym= x(:,k);
    % the delayed tail of prevSym will be superposed on the CP of currSym
    % resulting in inter-symbol interference
    concat = [prevSym.',currSym.'];
    temp=sum((current_rayleigh2.*delay(concat.',d,fs))*g',2);
    y(:,k) = temp(InRows+1:end);
    impulse_res(:,k) = sum((current_rayleigh1.*delay(delta,d,fs))*g',2);
end

% Downsample the output of the channel to match the input sample rate.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function used to decimate channel samples by a factor of two.

function y=ofdm_resample(x,type)

% Only need to calculate the impulse response once. This filter applies to
% all FFT size options. Specs are derived from the worst case (2048 pt).

persistent h
if isempty(h)
    Nused=852;
    Nfft=2048;
    h=remez(51,[0 (Nused+1)/Nfft 1-(Nused+1)/Nfft 1],[1 1 0 0],[1 150]);
end

% Use the resample function since it removes the delay.

if type == 'up'
    y=resample(x,2,1,h);
elseif type == 'dn'
    y=resample(x,1,2,h);
end