function OutSig = SimulateCOSTChannel( InSig, ChannelParams, fs)
% SimulateCOSTChannel - simulate multi-path fading channel based on COST
%                       259 model
%
% Compute the output signal from a multi-path fading channel from a given
% input signal and the model for the channel.
%
% Inputs:
%   InSig          - baseband equivalent input signal
%   ChannelParams  - structure ChannelParams must have fields
%                   Delay   - relative delay
%                   Power   - relative power in dB
%                   Doppler - type of Dopller spectrum
%                   fd      - max. Doppler shift
%   fs             - sampling rate
%
% Outputs:
%   OutSig        - baseband equivalent output signal
%
% Example:
%   Channel    = tux;
%   Channel.fd = 30;
%   OutSig = SimulateCOSTChannel( InSig, ChannelParams, fs)

%% convert powers to linear scale
Power_lin = dB2lin( ChannelParams.Power);

%% Bin the delays according to the sample rate
QDelay = floor( ChannelParams.Delay*fs );

% set surrogate delay for each bin, then sum up the power in each bin
Delays = ( ( 0:QDelay(end) ) + 0.5 ) / fs;
Powers = zeros( size(Delays) );
for kk = 1:length(Delays)
    Powers( kk ) = sum( Power_lin( QDelay == kk-1 ) );
end

%% Time-varying channel coefficients
% the channel taps are slowly time-varying. We keep them constant until the
% phase change for a tap may have reached phi degrees based on maximum
% Doppler shift.
maxPhi  = 0.1;     % maximum phase shift before filter taps are updated
if (ChannelParams.fd == 0)
    DeltaH = length(InSig)/fs;
else
    DeltaH  = maxPhi/360 / ChannelParams.fd; % time until maxPhi phase change
end

QDeltaH = round(DeltaH*fs);              % same in sample periods

% guard against very small sampling rates
if (QDeltaH < 2)
    QDeltaH = 2;
    DeltaH = QDeltaH/fs;
end

% DeltaH is the sampling rate of the filter used to update channel
% coefficients. The frequency response of that filter is the (square root
% of) the Doppler spectrum.
% DeltaH is generally to small to give us meaniningful frequency
% resolution; we downsample here and upsample later. The downsampling
% factor provides a sampling rate that is approximately 10 times the
% Doppler frequency
if (ChannelParams.fd == 0) % special case
    Down = 1;
    hh = 1;
else
    Down   = ceil( 1/(10*DeltaH*ChannelParams.fd) );
    NH     = 256; % number of taps for filter
    ff = (0:NH-1)/((NH-1)*Down*DeltaH);

    % Frequency response and impulse response of the filter
    % Fix Me: support other Doppler types
    % desired frequency response of filter:
    HH   = sqrt( ClassDoppler( ff, ChannelParams.fd ) );
    % design filter with desired frequency response
    hh = Persistent_firpm( NH-1, 0:1/(NH-1):1, HH );
    hh = hh/norm(hh);   % ensure filter has unit norm
end

% Compute enough filter taps to cover the signal
NumPeriods = ceil( length(InSig)/QDeltaH );

% array to keep filter taps
Taps = zeros( NumPeriods, length(Powers) );

% we need NumPeriods samples for each tap, we ensure that we are past the
% transient at the beginning
NSamples = length( hh ) + max( ceil(NumPeriods/Down), 10 );

for kk = 1:length(Powers)
    % generate a white Gaussian random process
    ww = sqrt( Powers( kk )/2)*...
        ( randn( 1, NSamples) + j*randn( 1, NSamples) );
    % filter so that spectrum equals Doppler spectrum
    ww = conv( ww, hh );
    ww = ww( length( hh )+1:NSamples ).';
    % interpolate to a higher sampling rate 
    ww = interp( ww, Down );
    % ww = interpft(ww, Down*length(ww));
    % store time-varying filter taps for later use
    Taps( :, kk ) = ww(1:NumPeriods);
end

%% Filter input signal
% we call conv repeatedly, each time with different filter coefficients
% as stored in Taps. 

% storage for output signal
OutSig = zeros(1, length(InSig)+length(Powers)-1);
% initialize loop variables
Start = 1;
kk = 1;

while ( Start < length(InSig) )
    EndIn  = min( Start+QDeltaH, length(InSig) );
    EndOut = EndIn + length(Powers)-1;
    OutSig(Start:EndOut) = OutSig(Start:EndOut) + ...
        conv( Taps(kk,:), InSig(Start:EndIn) );
    
    kk = kk+1;
    Start = EndIn+1;
end



%% Private functions
function HH = ClassDoppler(ff, fd)
% ClassDoppler - classical Doppler spectrum

HH = zeros( size( ff ) );

isSmall = ( abs(ff) < fd );

HH( isSmall ) = 1./sqrt( 1 - (ff(isSmall)/fd).^2 );



function hh = Persistent_firpm( order, FF, HH )
% Persistent_firpm - avoid computing the same filter all the time

persistent known_order known_FF known_HH known_hh

if ( all(size(known_order) == size(order) ) && known_order == order && ...
        all(size(known_FF) == size( FF)) && all( known_FF == FF ) && ...
        all(size(known_HH) == size( HH)) && all( known_HH == HH ) )
    hh = known_hh;
else
    hh = firpm( order, FF, HH );
    known_order = order;
    known_FF    = FF;
    known_HH    = HH;
    known_hh    = hh;
end
