How can I implement a stable Shelving filter for short pulses?

I don't know much about MATLAB and have only been experimenting with it for 1 week.

I have a SOFA file with a HRTF (Head Related Transfer Function) measurement.

The measurement is faulty and I wanted to apply a simple shelf filter below 350Hz on the left channel. After a lot of research I managed to apply an IIR shelf filter with the "shelvFilt" function. The problem is that the frequency response is very unstable, probably because HRTFs are very short pulses.

The next idea that comes to my mind is an FIR filter, does anyone have an idea how to best implement this?

This is the code with the IIR shelf filter using the SOFAToolbox + a screenshot with the unstable result.

% Load the SOFA file
data = SOFAload('/Users/username/Desktop/HRTF_48000.sofa');
% Extract the left channel data from Data.IR
leftChannelData = squeeze(data.Data.IR(:, 1, :));
% Create the shelving filter
shelvFilt = shelvingFilter(Gain=-1, Slope=1, CutoffFrequency=350, FilterType="lowpass", SampleRate=48000);
% Apply the shelving filter to the left channel data
filteredLeftChannelData = shelvFilt(leftChannelData);
% Update the left channel data in the new SOFA structure
newData = data;
newData.Data.IR(:, 1, :) = filteredLeftChannelData;
% Save the new SOFA structure as a new SOFA file
SOFAsave('/Users/username/Desktop/HRTF_48000_Modified.sofa', newData, 0);

[Paul]

Hi Hasan,

Is this the desired magnitude response fo the filter you want to apply to the input signal?

shelvFilt = shelvingFilter(Gain=-1, Slope=1, CutoffFrequency=350, FilterType="lowpass", SampleRate=48000);
[b,a] = shelvFilt.coeffs;
[h,f] = freqz(b,a,1024,48000);
semilogx(f,db(abs(h))),grid
ylabel('Mag (db)')

Can you expand on what you mean by " frequency response is very unstable"?

[Hasan Aydin]

Hi Paul,

Thanks!

That's exactly how it should look like, maybe with a little faster slope, so that the signal is attanuated at 200Hz.

the unstable frequency response you can see in my first screenshot i attached in my post, but it should look like a clean shelving like the plot you sent me.

You can also take a look at the screenshot I am attaching now, maybe this helps

[Paul]

How did you generate the frequency repsponse plot attached to the question?

How are you generating the lower graph in the plot attached to the preceding comment? Using the coded in the question?

Can you attach the left and right channel data in .mat file using the paperclip icon in the Insert menu?

[Hasan Aydin]

I created the diagrams with audio analysis software like REW and Plugin Doctor.

In the first screenshot, I calculated the difference between the original SOFA file response and the matlab processed one in the REW software.

In the second screenshot, I used a frequency analyser for the bottom plot, which analyses the output of a spatial renderer software. It renders the SOFA file containing the pulses, I think over 1500 pulses for different directions/head positions.

My code basically extracts all the pulses (from the left channel), processes them all at once and puts them back into a SOFA container.

This is where I'm having trouble, because I wanted to process all the impulses with a clean low shelf filter. But after importing the processed SOFA file into Spatial Renderer and calculating the difference from the original SOFA file in the REW software, I get this weird looking low shelf plot...

Since I have almost no skills with Matlab, this is the only way I can do it. But I can send you the SOFA file if you are familiar with the SOFAToolbox in Matlab...

[Paul]

I have no idea what the SOFA toolbox is, nor the REW software, nor the Spatial Renderer. As this is a Matlab forum, I think the best you can hope for is to save leftChannelData and rightChannelData to a .mat file and attach the .mat file to the Question or a comment using the paperclip icon so that people can take a look at it based on the code included in the Question.

[jibrahim]

Hi Hasan,

If you attach your SOFA file, I can try to take a look.

Two comments:

1) Starting in release R2023B, you will be able to read and write SOFA files in MATLAB with Audio Toolbox.

2) Make sure leftChannelData is a column vector, as that is what the shevling filter expects. A row vector would be wrong, as each sample of the input would be interpreted as coming from an independent channel.

[Hasan Aydin]

Hi Ibrahim,

thank you for this hint!

I have now tried to make a column vector (with the help of ChatGPT). The code is executed without errors but still looks like in my first screenshot.

I am uploading my SOFA file, maybe you can look over it and help me.

Unfortunately the file is 9.5MB, too big to upload directly, but I created a download link for you on WeTransfer, I hope that's ok for you?...

https://we.tl/t-lhX8F0YFq2

Thanks a lot!

% Load the SOFA file
data = SOFAload('/Users/hasanaydin/Desktop/HRTF_DF_Hasan_ITD_encoded_48000.sofa');

% Extract the left channel data from Data.IR
leftChannelData = squeeze(data.Data.IR(:, 1, :));

% Create a shelving filter
shelvFilt = shelvingFilter(-1.1, 1.5, 350, "lowpass");

% Initialize the filtered data array
filteredData = zeros(size(leftChannelData));

% Apply the shelving filter to each time slice of the left channel data
for i = 1:size(leftChannelData, 3)
    filteredData(:, :, i) = shelvFilt(leftChannelData(:, :, i));
end

% Update the left channel data in the new SOFA structure
newData = data;
newData.Data.IR(:, 1, :) = filteredData;

% Save the new SOFA structure as a new SOFA file
SOFAsave('/Users/hasanaydin/Desktop/HRTF_DF_Hasan_ITD_encoded_48000_modified.sofa', newData, 0);

Answer 1

[jibrahim]

Hi Hasan,

Since you are trying to modify a short impulse response, you might want to apply the transformation in the frequency domain, and then go back to the time domain. The following code takes advantage of R2023B functionality (reading SOFA files). The basic idea should work in any release though:

% Design your shelving filter
shelvFilt = shelvingFilter(Gain=-5, Slope=1,...
    CutoffFrequency=350, FilterType="lowpass", SampleRate=48000);

% View the filter response
visualize(shelvFilt)

% You can change filter parameters and observe the change in the response
shelvFilt.Gain=-1;

% Get the frequency response of the shelving filter
[b,a] = coeffs(shelvFilt);
NFFT = 1024;
HShelv = freqz(b,a,NFFT);

% Read the SOFA data
s = sofaread('HRTF_DF_Hasan_ITD_encoded_48000.sofa');

% Get one impulse response
imp = squeeze(s.Numerator(10,1,:));
impLen = length(imp);
H = fft(imp,NFFT);

% Apply shelving filter
HNew = H.*HShelv;

% Back to time domain
impNew = ifft(HNew);
impNew = real(impNew(1:impLen,:));

figure
subplot(211)
plot(imp)
subplot(212)
plot(impNew)

figure
freqz(imp,1,NFFT)
hold on
freqz(impNew,1,NFFT)

Once you are happy with your filter, you cabn apply the change to all your impulse responses and generate a new SOFA file. For example:

% Read the SOFA data
s = sofaread('HRTF_DF_Hasan_ITD_encoded_48000.sofa');

% Design your shleving filter
shelvFilt = shelvingFilter(Gain=-5, Slope=1, CutoffFrequency=350, FilterType="lowpass", SampleRate=48000);
visualize(shelvFilt)

% Get the frequency response of the shelving filter
[b,a] = coeffs(shelvFilt);
NFFT = 1024;
HShelv = freqz(b,a,NFFT);

numMeas = size(s.Numerator,1);
impLen = size(s.Numerator,3);
for index=1:numMeas

    % Get an impulse response for left channel
    imp = squeeze(s.Numerator(index,1,:));
    H = fft(imp,NFFT);

    % apply shelving filter
    H = H.*HShelv;

    impNew = ifft(H);
    impNew = real(impNew(1:impLen,:));

    s.Numerator(index,1,:) = impNew.';
end

% Write new SOFA file
sofawrite('NewData.sofa',s);

If you want to apply a filter on top of your original HRTF data in a streaming loop, you can do something like this:

% Read the SOFA data
s = sofaread('HRTF_DF_Hasan_ITD_encoded_48000.sofa');

% choose your impulse response
imp = squeeze(s.Numerator(1,1,:));
fir = dsp.FrequencyDomainFIRFilter(imp.');

% Design your shleving filter
shelvFilt = shelvingFilter(Gain=-5, Slope=1, CutoffFrequency=350,...
                 FilterType="lowpass", SampleRate=48000);
visualize(shelvFilt)
% Use this UI to tune the filter as the sim is running
parameterTuner(shelvFilt)

afr = dsp.AudioFileReader('FunkyDrums-48-stereo-25secs.mp3');
adw = audioDeviceWriter(SampleRate=48000);
spect = spectrumAnalyzer(SampleRate=48000,PlotAsTwoSidedSpectrum=false,FrequencyScale="log");
while ~isDone(afr)
    frame = afr();
    x = fir(frame);
    y = shelvFilt(x);
    adw(y);
    spect(y(:,1))
    drawnow limitrate
end

[Hasan Aydin]

Hi Ibrahim,

you helped me a lot with this! your code didn't work right away but I adapted my old code accordingly and experimented, also learned new things. The key was instead of applying the shelvFilt directly, to apply the coefficients of it. That worked exactly the way I wanted it to!

Thanks a lot!

% Load the SOFA file
data = SOFAload('/Users/hasanaydin/Desktop/HRTF_DF_Hasan_ITD_encoded_48000.sofa');

% Extract the left channel data from Data.IR
leftChannelData = squeeze(data.Data.IR(:, 1, :));

% Create a shelving filter
shelvFilt = shelvingFilter(-1, 1, 350, "lowpass");

% Get the filter coefficients
[b, a] = coeffs(shelvFilt);

% Apply the shelving filter to the left channel data
filteredLeftChannelData = filter(b, a, leftChannelData, [], 2);

% Update the left channel data in the new SOFA structure
newData = data;
newData.Data.IR(:, 1, :) = filteredLeftChannelData;

% Save the new SOFA structure as a new SOFA file
SOFAsave('/Users/hasanaydin/Desktop/HRTF_DF_Hasan_ITD_encoded_48000_modified.sofa', newData, 1);

[Paul]

Is it not odd that shelvFilt couldn't be applied directly? It seems that's supposed to be basic functionality of the shelvingFilter.

[jibrahim]

You should be able to apply it directly:

 s = shelvingFilter(Gain=-1);
[b,a] = coeffs(s);
in = randn(1024,1);
y1 = s(in);
y2 = filter(b,a,in);
isequal(y1,y2) % true

[Hasan Aydin]

but here you also get the coefficients of this filter with [b,a] = coeffs(s); or am I wrong? What I meant by using shelvFilt directly is as in my first code, without using [b,a] = coeffs at all.

Maybe I didn't fully understand what coeffs of a filter are :)

[jibrahim]

Hi Hasan, you are getting the filter coefficients correctly. I am just pointing out that using the coefficients with the filter function should be exatcly the same as using the objhect directly.