1
$\begingroup$

I am trying to make a frequency domain adaptive filter in matlab. It uses matlab adaptfilt.fdaf to create the filter parameters like step size and initializing initial filter weight values. Then I have tried to implement the overlap - save frequency domain adaaptive filter algorithm from the paper "Frequency-domain and multirate adaptive filtering" by J. J. Shynk. It is also based on the source code of the thisfilter.m script that can be found in "%MATLAB_FOLDER%\R2009b\toolbox\filterdesign\flterdesign\@adaptfilt\@fdaf\thisfilter.m". I used sine wave as the desired signal and

noisy = sinewave.*(rand(1, length(t) ).') * 0.1;  

as the noisy input. At the end, I plot the difference (output - desired) to see how close the filter output comes to the desired signal. For the above noisy input, if I play around witht he step size and other filter parameters, I can make this difference go closer to 0 (although it never actually becomes 0, but sort of oscillates around 0) as I process more and more blocks of the input signal. However, if I then change the input to

noisy = sinewave + (rand(1, length(t) ).') * 0.1;

I just can't make the difference converge towards 0 no matter what I do. Either the difference increases to very high magnitude (around the range x$10^{50}$) if the step size is too high, or the difference oscillates in a roughly sine wave pattern, but never seems to go towards 0. Can anyone provide any ideas on how this can be fixed?

A few facts that I should point out here:
When the noisy signal is noisy = sinewave.*(rand(1, length(t) ).') * 0.1;, the plot of the difference (output - desired) is in the range +0.3 to -0.3 for a good value of step size, that is, one that makes this plot converge and not diverge. When the noisy signal is noisy = sinewave + (rand(1, length(t) ).') * 0.1; this difference is of the range +4 to -4 when the output does not diverge.

If I get some good suggestions on how I can improve this code I will post further details about the source code and output.
Here is the current version of the Matlab code. It still has some lines of code that I used for my own trouble shooting so it's not exactly like Shynk's paper right now, but this version can be compiled. It uses the noisy = sinewave.*(rand(1, length(t) ).') * 0.1; version of the noisy signal. If I get some good suggestions from here I will update the source code as well for further discussion.

clc
clear all
close all
% % %  Arguments to fdaf % % %
% H = ADAPTFILT.FDAF(L,STEP,LEAKAGE,DELTA,LAMBDA,BLOCKLEN,OFFSET,COEFFS, STATES)
L = 128;            % number of filter coefficients ( = n for best results; N = Block length; default = 10)
% L is the number of samples by which shifting takes place during every
% iteration
STEP = 0.6;    % mu; default = 1; corresponds to 2*mu in paper
LEAKAGE = 1;    % filter leakage, default = 1 --> No leakage
DELTA = 1;      % initialize and assign FFT input signal powers; should be 1 if it's effect is to be ignored; default = 1
LAMBDA = 1;     % assign averaging factor; should be 1 if it's effect is to be ignored; default = 1
BLOCKLEN = 128; % length of block, number of samples read from the input to be filtered at one iteration --> N; should equal L for best efficiency
OFFSET = 0;     % should be 0 to remove it's effect; default = 0
COEFFS = zeros(L, 1);   % probably initial value for filter coefficients; should be 0; default = 0 array of length L
STATES = zeros(L, 1);   % the very first "previous" value for input noisy signal

% % % % % % % % % % Define input signal % % % % % % % % % % 
%% Time specifications:
Fs = 8000;                   % samples per second
dt = 1/Fs;                   % seconds per sample
StopTime = 1.6;             % seconds
t = (0:dt:StopTime-dt)';     % seconds
%% Sine wave:
Fc = 60;                     % hertz
sinewave = cos(2*pi*Fc*t);
desired = sinewave;
% load rand_seq.dat
rand_seq = (rand(1, length(t) ).') * 0.1;
noisy = sinewave.*rand_seq;   % noisy signal

% % % % % % % % % % Creating filter % % % % % % % % % % 
h = adaptfilt.fdaf(L,STEP,LEAKAGE,DELTA,LAMBDA,BLOCKLEN,OFFSET,COEFFS, STATES);
x = noisy;
d = desired;
% % % % % % % % % % Applying filter % % % % % % % % % % 
ntr = length(x);            %  temporary number of iterations 
N = h.BlockLength;          %  block length
L = h.FilterLength;         %  number of coefficients
ntrB = floor(length(x)/N);  %  temporary number of block iterations ; --> has to be exactly divisible in the original function
y = zeros(size(x));         %  initialize output signal vector
e = y;                      %  initialize error signal vector
X = zeros(L+N,1);           %  initialize temporary input signal buffer
E = zeros(L+N,1);           %  initialize temporary error signal buffer
Ef = zeros(L+N,1);          %  initialize temporary error signal buffer
nnL = 1:L;                  %  index variable used for input signal buffer
nnLpN = N+1:N+L;            %  index variable used for input signal buffer
nnNpL = L+1:L+N;            %  index variable used for input signal buffer
nnLpNr = L+N:-1:N+1;        %  index variable used for coefficient updates
FFTW = h.FFTCoefficients.'; %  initialize and assign frequency domain Coefficients
normFFTX = h.Power;         %  initialize and assign FFT input signal powers
X(nnLpNr) = h.FFTStates;    %  assign input signal buffer
mu = h.StepSize;            %  assign step size
ZN = zeros(N,1);            %  assign N-dimensional zero vector
mu_k = zeros(L+N, 1);       % use seperate mu value for each frequency bin in the fft domain
                            % mu gets multiplied to values in the FFT
                            % domain, but each element of mu is a scalar
Pest = zeros(L+N, 1);
lambda = 0.8;              
alpha = 1 - lambda;

FFTW_sample_index = zeros(ntrB, 1);

for n=1:ntrB,
    nn = ((n-1)*N+1):(n*N); %  index for current signal blocks
    X(nnL) = X(nnLpN);      %  shift temporary input signal buffer up
    X(nnNpL) = x(nn);       %  assign current input signal vector
    FFTX = fft(X);          %  compute FFT of input signal vector; FFTX is the result of an L+N point FFT

%     if n == 1
%         Pest = (abs(FFTX)).^2;
%     else
%         Pest = lambda * Pest + alpha * ( (abs(FFTX)).^2 );
%     end

%     mu_k = mu ./ Pest;

%     if mod(n, 30) == 0
% %         mu = 0.009;
%         FFTW(:, 1) = 0;
%     end

    mu_k = mu ./( lambda + (abs(FFTX)).^2 );
    Y = ifft(FFTW.*FFTX);   %  compute current output signal vector
    y(nn) = Y(nnNpL);       %  assign current output signal block
    e(nn) = d(nn) - y(nn);  %  assign current error signal block
%     E(nnNpL) = mu*e(nn);    %  assign current error signal vector
    E(nnNpL) = e(nn);
    FFTE = fft(E);          %  compute FFT of error signal vector
%     normFFTX = bet*normFFTX + ombet*real(FFTX.*conj(FFTX));  
                            %  update FFT input signal powers
%     G = ifft(FFTE.*conj(FFTX)./(normFFTX + Offset));  
                            %  compute gradient vector

    FFTE_multConjFFTX = FFTE.*conj(FFTX);
    G = mu_k .* FFTE_multConjFFTX;
    G = ifft(G);

%     G = ifft(FFTE.*conj(FFTX));
    G(nnNpL) = ZN;          %  impose gradient constraint
%     FFTW = lam*FFTW + mu * fft(G);  %  update frequency domain coefficients
%     FFTW = FFTW + mu * fft(G);    % modified to FFTW = FFTW + mu_k .* fft(G);
%       FFTW = FFTW + mu_k .* fft(G);
    FFTW = FFTW + fft(G);       % USe this equation if mu_k .* FFTE_multConjFFTX; has been used above

    FFTW_sample = abs(FFTW(25));
%     invFFTW_index(n) = sum( invFFTW )/ntrB;
    FFTW_sample_index(n) = FFTW_sample;
% 
%     if n < 20
%         FFTW = FFTW + fft(G);       % USe this equation if mu_k .* FFTE_multConjFFTX; has been used above
%     end

end

if isreal(x) && isreal(d),
    y = real(y);            %  constrain output signal to be real-valued
    e = real(e);            %  constrain error signal to be real-valued
end

diff_here = y - desired;
figure(1); plot(diff_here);
figure(2); plot(1:ntrB, FFTW_sample_index);
% figure(2);
% plot(1:length(noisy), noisy, 'b', 1:length(y), y/4, 'r');
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.