%% Estimate Coefficients of long FIR Filter using Frequency-Domain Adaptive Filter
% Use a frequency-domain adaptive filter to estimate the coefficients of a long 
% FIR filter. The FIR filter models the impulse response of a room. Use the partitioned 
% mode in the frequency-domain adaptive filter to reduce filter latency.
% 
% *Note:* This example runs only in R2018a or later.
%% Initialization
% Generate a long FIR impulse response of 8192 samples and assign the impulse 
% response to a |dsp.FIRFilter| object, |room|. This object models the impulse 
% response of a room. Create a |dsp.FrequencyDomainAdaptiveFilter| filter, |lmsfilt|, 
% in partitioned constrained mode. Set the length of the filter to one-fourth 
% the length of the impulse response of the room. Set the block length of the 
% filter to 128 samples. Set the step size to 0.025, initial power to 0.01, averaging 
% factor to 0.98, offset to 1, and the leakage factor to 1. Initialize a |dsp.ArrayPlot| 
% object to view the filter coefficients. Initialize a |dsp.TimeScope| object 
% to show the mean-squared error between the filter output and the desired signal.
%%
fs = 16e3;
m = 8192;
[b,a] = cheby2(4,20,[0.1 0.7]);
impulseResponseGenerator = dsp.IIRFilter(...
    'Numerator', [zeros(1,6) b], ...
    'Denominator', a);
roomImpulseResponse = impulseResponseGenerator( ...
    (log(0.99*rand(1,m)+0.01).*sign(randn(1,m)).*exp(-0.002*(1:m)))');
roomImpulseResponse = roomImpulseResponse/norm(roomImpulseResponse);
room = dsp.FIRFilter('Numerator',roomImpulseResponse');

lmsfilt = dsp.FrequencyDomainAdaptiveFilter(...
    'Method','Partitioned constrained FDAF',...
    'Length',m/4, ...
    'BlockLength',128,...
    'StepSize',0.025, ...
    'InitialPower',0.01, ...
    'AveragingFactor',0.98,...
    'Offset',1,...
    'LeakageFactor',1);

FrameSize = lmsfilt.BlockLength; NIter = 2000;

AP = dsp.ArrayPlot('YLimits',[-0.2 .2],'ShowLegend',true, ...
    'Position',[0 0 560 420],'ChannelNames', ...
    {'Actual Coefficients','Estimated Coefficients'});

TS = dsp.TimeScope('SampleRate',fs,'TimeSpan',FrameSize*NIter/fs,...
    'TimeUnits','Seconds',...
    'YLimits',[-50 0],'Title','Learning curve',...
    'YLabel','dB', ...
    'BufferLength',FrameSize*NIter,...
    'ShowGrid',true);

signalmean = dsp.MovingAverage('SpecifyWindowLength',false);
                         
%% Streaming
% Generate a random input signal using the |randn| function. The frame size 
% of the input matches the block length of the adaptive filter. The desired signal 
% is the sum of the output of the FIR filter (room) and a white Gaussian noise 
% signal. Pass the input signal and the desired signal to the adaptive filter. 
% Compute the adaptive filter output and the error between the output and the 
% desired signal.
% 
% Estimate the time-domain coefficients of the adaptive filter by taking 
% the IFFT of the frequency-domain coefficients vector returned by the |lmsfilt.FFTCoefficients| 
% property. Compare the estimated coefficients with the actual coefficients assigned 
% to the FIR filter (room). Once the adaptive filter has converged its output 
% to the desired signal, and minimized the error signal, the estimated coefficients 
% match closely with the actual coefficients. This means that the adaptive filter 
% has successfully adapted itself to model the impulse response of the FIR filter 
% (room).
%%
for k = 1:NIter
    x = randn(FrameSize,1);                
    d = room(x) + 0.01*randn(FrameSize,1); 
    [y,e] = lmsfilt(x,d);
    FFTCoeffs = lmsfilt.FFTCoefficients;
    w = ifft(FFTCoeffs,[],2,'symmetric');
    w = w(:,1:FrameSize) + w(:,FrameSize+1:end);
    w = reshape(w.',1,m/4);
    AP([roomImpulseResponse(1:m/4),w.']);
    TS(10*log10(signalmean(abs(e).^2)));
end
%%
% 
% <<../arrayplot_adaptivefilter.png>>
% 

%% 
% As the filter adapts with time, you can see in the time scope that the 
% mean-squared error becomes minimal. Simultaneously, the estimated coefficients 
% match the actual coefficients closely in the array plot.

%% 
% Copyright 2012 The MathWorks, Inc.