With help from the dsp guide and DSP-Related, I tried to implement a leaky integrator in octave/matlab. It seems to work in general, but there are a few problems.

So that's what I do:

  1. Calculate the optimal integrator (1/jω, damping below fc Hz) in the frequency domain
  2. Band limit this filter (using custom_filter.m)
  3. Apply the overlap-add method

The result can be seen in the images Fig 1 to Fig 4. As Fig 3 shows, the band limiter makes a linear phase shift from the desired constant phase shift, due to shift in time domain (see custom_filter.m). So my integrator isn't working, because the phase shift is wrong.

There's an explanation in the update below the code

No distortion, wrong phase response

If I undo the shift in time domain, I get a filter that looks quite good (phase shift is at least constant...), but the resulting signal looks very distorted (second image).

Distortion, but better phase response

So where does this distortion come from? And more important, how can I avoid it?

So here's the code. I was hoping I could attach it as a zip but couldn't find how...


% --------------------------------------------------------------------------------
% Helper functions, used below. All in one file to copy and paste it
% to stackexchange.

#  y = spec_integrate( f, fc )
#  Integrates the f Hz line of a complex spectrum
#  with a fc Hz highpass
function y=spec_integrator(f, fc)
  if f==0,
    y=1* 1000/( i * 2*pi*f );

    # attenuate at f < fc
    if f < fc,
    y = y*d;

% y = ccmirror( x ) "complex conjugate mirror" where:
% length( x ) = 2^N
% x contains the positive frequencies of a spectrum
% x(1    )   is DC
% x(N/2+1)   is the folding frequency
% x((N/2+2):N) is ignored
% y is created from x so that its negative frequencies y((N/2+2):N)
% are the complex conjugate of the positive frequencies.
% This means imag(ifft( y ) ) = 0
function y = ccmirror( x )
   y = x;
   y((N/2+2):N) = conj(y((N/2):-1:2));

% [impulse_response actual_frequency_response] = custom_filter( desired_frequency_response, M)
% Band limit a filter that has been customized in frequency domain
% desired_frequency_response is the custom filter definition
% M is the size of the filter kernel / impulse response
% impulse_response is the bandlimited impulse response, for
%                  informational purpose.
% actual_frequency_response is the bandlimited frequency response
% This is taken from [DSP-Related][2]
function [impulse_response actual_frequency_response] = custom_filter( desired_frequency_response, M)

  % Get the aliased impulse response
  impulse_response=ifft( desired_frequency_response );
  % Rotate the signal M/2 points to the right
  % Note: translation in time domain means a phase shift in frequency
  % domain! The phase of the frequency response is modified!
  impulse_response=circshift( impulse_response, [0,M/2] );

  % Truncate and window the signal
  impulse_response=impulse_response .* [ hanning(M)' zeros(1,N-M)];

  %% One could undo the phase shift from above ... but it gives ugly
  %% results. Why?
  %% impulse_response=circshift( impulse_response, [0,-M/2] );


% --------------------------------------------------------------------------------
%  Now comes the more interesting stuff
% --------------------------------------------------------------------------------

% Set up simulation parameters

L = 1024;                                 % FIR filter length in taps
fc = 10;                                  % lowpass cutoff frequency in Hz
fs = 4000;                                % sampling rate in Hz

Nsig = 8192;                               % signal length in samples
fsig = 159;                               % Frequency of the demo signal
asig = 1;                                 % Amplitude of the demo signal

%FFT processing parameters:

M = L;                                    % nominal window length
Nfft = 2^(ceil(log2(M+L-1)));             % FFT Length
M = Nfft-L+1;                             % efficient window length
R = M;                                    % hop size for rectangular window
Nframes = 1+floor((Nsig-M)/R);            % no. complete frames

lc = round(fc * Nfft/fs);                 % line of the corner frequency

% Generate the test signal:

sig = asig * sin( 2 * pi * (0:(Nsig-1)) * fsig/fs );

% Define the filter in frequency domain

% integrator
H = zeros( 1, Nfft );
for f = 1:(Nfft/2+1)
  H(f) = spec_integrator( (f-1) *fs/Nfft, fc );
H = ccmirror( H );

Hoptimal = H;

[impulse_response Hactual] = custom_filter( Hoptimal, L);
H = Hactual;

% Carry out the overlap-add FFT processing:

y = zeros(1,Nsig + Nfft);                 % allocate output+'ringing' vector
for m = 0:(Nframes-1)
    index = m*R+1:min(m*R+M,Nsig);        % indices for the mth frame
    xm = sig(index);                      % windowed mth frame (rectangular window)
    xmzp = [xm zeros(1,Nfft-length(xm))]; % zero pad the signal
    Xm = fft(xmzp);
    Ym = Xm .* H;                         % freq domain multiplication
    ym = real(ifft(Ym));                  % inverse transform
    outindex = m*R+1:(m*R+Nfft);
    y(outindex) = y(outindex) + ym;       % overlap add

f=(0:(Nfft-1))*fs/Nfft;                   % frequencies for plotting
tin=(0:(length(sig)-1))/fs;               % times for plotting

figure 1 % input
plot( tin, sig );
hold on;
for m = 1:Nframes
  plot( [ m*L m*L ]/fs, [ min(sig), max(sig) ], "r");
title( "Fig 1: Input signal")
legend( sprintf("%g Hz", fsig), "Overlap-Add-Frames");

figure 2 % transfer function (abs)
loglog( f, abs(H), "r" );
hold on;
loglog( f, abs(Hoptimal));
hold off;
title( "Fig 2: Filter magnitude")
legend( "Band limited", "Optimal");

figure 3 % transfer function (phase)
plot( f, unwrap(arg(H)), "r" );
hold on;
plot( f, unwrap(arg(Hoptimal)));
hold off;
title( "Fig 3: Filter phase")
legend( "Band limited", "Optimal");


figure 4 % result
plot( tout, y );
hold on;
plot( [ L/2 L/2 ]/fs, [ min(y) max(y) ], "r");
plot( [ L L ]/fs, [ min(y) max(y) ], "r");
plot( [ length(y)-Nfft length(y)-Nfft ]/fs, [ min(y) max(y) ], "r" )
plot( [ length(y)-Nfft-L/2 length(y)-Nfft-L/2 ]/fs, [ min(y) max(y) ], "r" )
plot( [ length(y)-Nfft-L length(y)-Nfft-L ]/fs, [ min(y) max(y) ], "r" )
hold off;
title( "Fig 4: Filtered signal")
legend( sprintf("%g Hz", fsig), sprintf( "L/2=%d samples / %gs", L/2, L/2/fs));


Ok, so I was thinking a bit about the phase response in image 1, Figure 3: The bandlimiter shifts the signal in time domain, which means a change of the gradient of the phase response. That's what we see in Figure 3. I've been unsettled by this because I didn't know what this means for my integrator. But then I went on and thought about the group delay. The derivative of my phase response is constant, so there is a constant group delay which just means that the time domain signal is shifted by a constant value - which is just what the bandlimiter did :) So the integrator is in fact working and to be exact I'd have to shift the final signal back to where it belongs.

Great, so now I think the integrator works. Now only one thing is left: why is the result distorted when I undo the time shift of the bandlimiter?

  • $\begingroup$ Markus:"So the integrator is in fact working" - because the group delay is constant? Do you think this is a criterion of a working integrator? $\endgroup$ – LvW Sep 17 '14 at 7:50
  • $\begingroup$ Well, I think it works because magnitude and phase shift match the 1/jω - criterion, at least if I fix the group delay which seems doable because it's constant. I have to admit that much of my question describes something that's not directly related to my problem: why is the result distorted if I undo the time shift of the band limiter? $\endgroup$ – Markus Sep 17 '14 at 7:55

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.