dsl.m

ADSL经典仿真程序

%DSL Simulate a standard digital subscriber line.
% [Y, H, N] = DSL(X, channelName, P, S) returns the received signal
% Y (without channel noise), channel impulse response in H, and
% channel noise in N. The channel noise consists of NEXT noise 
% and AWGN of power P. 
%
% X is the transmitted signal. channelName is a string containing
% the channel name. The channel data should be in a subdirectory 
% "channel". Currently supported channels are the CSA loops from 1 
% to 8. The name format is "csaloop#" where # = {1,2,..,8}. P is the
% AWGN power. S is a subsampling factor. The default length of
% H is 512 samples. A shorter channel can be generated by subsampling
% H by a factor of S. 

% Copyright (c) 1999-2003 The University of Texas
% All Rights Reserved.
%  
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%  
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%  
% The GNU Public License is available in the file LICENSE, or you
% can write to the Free Software Foundation, Inc., 59 Temple Place -
% Suite 330, Boston, MA 02111-1307, USA, or you can find it on the
% World Wide Web at http://www.fsf.org.
%  
% Programmers:	Guner Arslan
% Version:        @(#)dsl.m	1.4	07/26/00
% 
% The authors are with the Department of Electrical and Computer
% Engineering, The University of Texas at Austin, Austin, TX.
% They can be reached at [e-mail].
% Guner Arslan is also with the Embedded Signal Processing
% Laboratory in the Dept. of ECE, http://signal.ece.utexas.edu.

function [y,h,noise,NEXTnoise,AWGN] = dsl(N,x,channelName,noisePower,bf,chbf)

% open a figure for progress bar

   [figHndl statusHndl] = setprogbar('Creating Channel Data ...');     

% set path for channel files 
addpath channels
% load channel data 
channelTime = ['load ',channelName,'.time'];
eval(channelTime)
% put channel data into h and time axis into t
string = ['t = ',channelName,'(:,1); h = ',channelName,'(:,2);'];
eval(string)
% sampling frequency
fsample = 2.208e6;
% normalization is required due to the way channel data is stored
h1 = h/2.208e6;
h1 = [h1 ;zeros(2048,1)];
% POTS splitter
[b a] = cheby1(5,1,4.8e3/fsample,'high'); %a=1; b =1;
% channel with splitter
h = filter(b,a,h1);
%h = decimate(h,5);
h = h(1:N);
%try to use a modified channel
if (chbf==1)
    %h=h-h(512);
    nStartTaper=floor(N*0.805); %%% floor(0.805*512)=412
    taperLength=N-nStartTaper;
    Taper=[ones(1,nStartTaper-1),[1:-1/taperLength:0]];
    h=h.*Taper';
end

%[b a] = butter(1,0.5);
%h = filter(b,a,h2);
%rampup = linspace(0,1,20);
%rampdown = linspace(1,0,20);
%h(1:20) = h(1:20).*rampup.';
%h(end-19:end) = h(end-19:end).*rampdown.';
%h(512:end) = h(512:end).*(length(h(512:end))-1:-1:0)'/(length(h(512:end)));
% channel with TEQ if TEQ is given
%hw = filter(wteq/(norm(wteq)),1,h2);
% filter transmit data with channel impulse response
%h = decimate(h,5);
%h = [h(1:425); zeros(512-425,1)];

filtX = filter(h,1,x);
% update progress bar

   updateprogbar(statusHndl,1,4);
 
% prepare a frequency axis for a filter prototype
MM = 128;
f = linspace(0,fsample/2,MM);
% desired frequency response of the NEXT coupling filter

users = 8;
fo = 276e3;
Kmask = zeros(size(f));
Kmask(f>=28e3 & f<138e3) = 113e3*0.1*10^(-38/10);
Kmask(f>138e3) = 113e3*0.1*10.^( (-38-24*(f(f>138e3) - 138e3)/43125) /10);
PSDadsl = Kmask.*sinc(f./fo).^2;
xn = 0.882*1e-14*users.^(0.6);
Hnext = PSDadsl.*xn.*f.^1.5;
Anext = 1;
%Hnext(Hnext<=eps) = eps;

% normalize to unity
f = f./max(f);  
% design IIR NEXT coupling filter with coeffs a and b
%[b a] = yulewalk(MM,f,sqrt(Hnext));
b = remez(128,f,sqrt(Hnext));
a = 1;

% update progress bar

   updateprogbar(statusHndl,2,4);

% filter a white noise signal with NEXT filter
% NEXTinput  = sqrt(cov(x))*randn(size(filtX));
NEXTinput = randn(size(filtX));
NEXTnoise =  filter( b,a, Anext*NEXTinput/sqrt(cov(NEXTinput)) );

% variance of the NEXT noise
covnext = cov(NEXTnoise);

% update progress bar

   updateprogbar(statusHndl,3,4);

% generate AWGN with power of AWGNpower
AWGN = sqrt(noisePower).*randn(size(filtX));
% variance of awgn
% covawgn = cov(awgn);
% add AWGN with NEXT noise
noise = NEXTnoise + AWGN; 
% noises passes through the equalizer
% noiseW = filter(wteq/(norm(wteq)),1,noise);
% add the noise to the channel filtered signal
y = filtX;
% update progress bar
%if bf == 1
   updateprogbar(statusHndl,4,4);
   % close progress bar
   close(figHndl);
%end