In real time filtering of a song, all the filter states seem to be NaN, so no sound plays. Not sure why...

First I compute the filter coefficients using this function

function [b, a] = computeCoeffs(Fs, Wp, Ws, Rp, Rs, ftype)
    % Normalize frequencies
    Wp = Wp / (Fs / 2);
    Ws = Ws / (Fs / 2);

    % Design the Butterworth filter order
    [n, Wn] = buttord(Wp, Ws, Rp, Rs);

    % Calculate filter coefficients
    [b, a] = butter(n, Wn, ftype);
end

And here it runs. I use frameCounter to ensure the while loop only plays the song for 10 seconds.

% Initialize Parameters
frameSize = 1024; % Buffer size

% Initialize audio devices
fileReader = dsp.AudioFileReader('IGottaFeeling.wav', 'SamplesPerFrame', frameSize, 'PlayCount', 1);
Fs = fileReader.SampleRate; % 16 kHz sample rate
deviceWriter = audioDeviceWriter('SampleRate', Fs); % Set up for mono output

% Set a duration for testing
durationInSeconds = 10; % Duration to play
totalFramesToPlay = floor(Fs * durationInSeconds / frameSize);
frameCounter = 0;

% Initialize gains in dB
gain1 = 0;
gain2 = 0;
gain3 = 0;
gain4 = 0;
gain5 = 0;

% Precompute filter coefficients for each band
[b1, a1] = computeCoeffs(Fs, 500, 800, 1, 50, 'low');
[b2, a2] = computeCoeffs(Fs, [500 2000], [300 2200], 1, 50, 'bandpass');
[b3, a3] = computeCoeffs(Fs, [2000 4000], [1800 4200], 1, 50, 'bandpass');
[b4, a4] = computeCoeffs(Fs, [4000 6000], [3800 6200], 1, 50, 'bandpass');
[b5, a5] = computeCoeffs(Fs, 6000, 5800, 1, 50, 'high');

% Initialize filter states for each channel
zi1 = []; 
zi2 = [];
zi3 = [];
zi4 = [];
zi5 = [];


% Create Spectrum Analyzer objects
saBefore = dsp.SpectrumAnalyzer('SampleRate', Fs, 'Title', 'Spectrum Before Filter');
saAfter = dsp.SpectrumAnalyzer('SampleRate', Fs, 'Title', 'Spectrum After Filter');

% Real-time Processing
while ~isDone(fileReader) && frameCounter < totalFramesToPlay
    % Read audio frame from file
    x = fileReader();
    
    % Convert stereo input to mono if necessary
    if size(x, 2) == 2
        x = mean(x, 2); % Convert stereo to mono by averaging the two channels
    end

    % Display the spectrum before filter
    %saBefore(x);
    
    % Convert dB gain to linear scale
    linearGain1 = 10^((gain1-0) / 20);
    linearGain2 = 10^((gain2+0) / 20);
    linearGain3 = 10^((gain3+0) / 20);
    linearGain4 = 10^((gain4+0) / 20);
    linearGain5 = 10^((gain5+0) / 20);     

    % Filter audio for each band and update filter states
    [ch1, zi1] = filter(b1, a1, x, zi1);
    [ch2, zi2] = filter(b2, a2, x, zi2);
    [ch3, zi3] = filter(b3, a3, x, zi3);
    [ch4, zi4] = filter(b4, a4, x, zi4);
    [ch5, zi5] = filter(b5, a5, x, zi5);
    
    % Combine channels     
    audioOut = ch1 ;%+ ch2 + ch3 + ch4 + ch5;

    % Display the spectrum after filter
    %saAfter(audioOut);

    % Output audio
    deviceWriter(audioOut);
    
    % Real-time gain adjustment code can go here
    frameCounter = frameCounter + 1;
end

release(fileReader);             % close the input file
release(deviceWriter); 

Answer 1

[Aishwarya]

Hi Samuel,

As per my understanding, you are performing filtering of an audio file in the above code, but the code produces NaN values in the filter state.

After reviewing the code, here are few ways that can help troubleshoot or resolve the issue:

• Filter Design Parameters: Consider that the filter design parameters (cutoff frequencies, passband, and stopband ripples) are within reasonable ranges and that the cutoff frequencies are less than the Nyquist frequency (half the sampling rate).
• Filter Order Too High: A high filter order can lead to coefficients that are sensitive to numerical errors. One suggestion is to change the values of “Rp” and “Rs” variables which are given as input to “buttord” function to get a lower order filter. Another suggestion is to add lower limit to the filter order in the “computeCoeffs” function definition. Here is an example code to perform the same:

    % Design the filter order
    [n, Wn] = buttord(Wp, Ws, Rp, Rs);
    
    % Ensure the filter order is not too high
    n = min(n, 6); % Limiting the order to 6 as an example
    
    % Calculate filter coefficients
    [b, a] = butter(n, Wn, ftype);
    

• Check Filter Coefficients: Ensure that the filter coefficients returned by “computeCoeffs” function are valid and do not cause the filter to become unstable. The stability of the filter can be checked by using “isstable” function inside the “computeCoeffs” function definition. Here is an updated version of “computerCoeffs” function:

    function [b, a] = computeCoeffs(Fs, Wp, Ws, Rp, Rs, ftype)
    % Normalize frequencies
    Wp = Wp / (Fs / 2);
    Ws = Ws / (Fs / 2);
   
    % Design the filter order
    [n, Wn] = buttord(Wp, Ws, Rp, Rs);
    
    % Ensure the filter order is not too high
    n = min(n, 6); % Limiting the order to 6 as an example
    
    % Calculate filter coefficients
    [b, a] = butter(n, Wn, ftype);
    
    % Check the stability of the filter
    assert(isstable(b, a), 'The designed filter is not stable.');
    end

• Analyze frequency response of the filter: Use the “fvtool” function to analyze the frequency response of the filter and check for stability and proper behaviour. Here is an example on how to use “fvtool” function to analyze the filter:

     fvtool(b2, a2);

If the frequency response looks unusual, try redesigning the filter with different specifications or consider different type of filter.

Please refer to the following MathWorks documentations for more information about the functions used:

• “buttord” function: https://www.mathworks.com/help/signal/ref/buttord.html
• “isstable” function: https://www.mathworks.com/help/signal/ref/isstable.html
• “fvtool” function: https://www.mathworks.com/help/signal/ref/fvtool.html

I hope this helps in resolving the issue.

Best Regards,

Aishwarya