Filter long audio data using cascading biquad filter

Question

I would like to lowpass filter a long pcm wav file. I designed a 20 order butterworth lowpass filter in Octave and transform the zero pole gain data to a series of second order sections. Then I am applying these SOS filters to each processing block on the long audio. The script is like this:

clear all; close all;
[y, fs, bit] = wavread('sweep.wav');

# design butter
[z, p, k] = butter(20, 1200 / (fs / 2), 'low');
[sos, k] = zp2sos(z, p, k);

# frame index
framesPerBuffer = 512;
currentFrame = 1;

# final condition buffer
zbuf(2) = 0;

while currentFrame + framesPerBuffer < size(y, 1)
    # read one frame
    data = y(currentFrame:currentFrame + framesPerBuffer - 1, 1);
    for i = 1:size(sos, 1)
        data = filter(sos(i, 1:3), sos(i, 4:6), data);
    end
    y(currentFrame:currentFrame + framesPerBuffer - 1, 1) = data;
    currentFrame += framesPerBuffer;
end
# apply gain to filtered signal
y = y * k;

The "sweep.wav" input is a mono pcm file with 48kHz sample rate and 16 bit depth with duration of 40 seconds. It is a step sweep signal with each sine tone duration about 1 second, frequency range from 200-12000 Hz. The overall wave shape is like this: After the filtering process above the output signal has the wavform of this: Since butterworth filter has the flat passband/stopband, the waveform is clearly not correct. When I listen to the filtered wave it is like a filtered wav plus some low frequency modulation.

I have also tried to add final condition z buffer for the filter function but it is still perform bad. Is there anything wrong for me filter process? How can I get the correct filtered data out of this long PCM? Thank you very much!

Answer 1

You are not storing the filter states between buffers.

In this case, each filter is a 2nd order filter, which has 2 state variables (the values stored in the 2 delays implementing each filter).

So for each new frame, you are starting with state 0, introducing a glitch between frames.

You must store the resulting state for each filter, and apply it for each next frame.

Built-in Function: [Y, SF] = filter (B, A, X, SI)

 If the fourth argument SI is provided, it is taken as the initial
 state of the system and the final state is returned as SF.  The
 state vector is a column vector whose length is equal to the
 length of the longest coefficient vector minus one.  If SI is not
 supplied, the initial state vector is set to all zeros.

So for each filter, you must store the corresponding SF, and provide as SI on the next frame.