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

Question

Hasan Aydin 2023-6-24

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1987718-how-can-i-implement-a-stable-shelving-filter-for-short-pulses

评论： jibrahim 2023-6-29

采纳的回答： jibrahim

Screenshot 2023-06-24 at 22.26.09.png

在 MATLAB Online 中打开

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);

7 个评论
显示 5更早的评论隐藏 5更早的评论

Hasan Aydin 2023-6-25

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...

jibrahim 2023-6-26

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 2023-6-26

编辑：Hasan Aydin 2023-6-26

在 MATLAB Online 中打开

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 2023-6-27

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1987718-how-can-i-implement-a-stable-shelving-filter-for-short-pulses#answer_1263583

在 MATLAB Online 中打开

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

5 个评论
显示 3更早的评论隐藏 3更早的评论

Hasan Aydin 2023-6-28

在 MATLAB Online 中打开

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);

Hasan Aydin 2023-6-28

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 2023-6-29

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.

请先登录，再进行评论。

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

7 个评论
显示 5更早的评论隐藏 5更早的评论

采纳的回答

5 个评论
显示 3更早的评论隐藏 3更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

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

7 个评论 显示 5更早的评论隐藏 5更早的评论

采纳的回答

5 个评论 显示 3更早的评论隐藏 3更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

7 个评论
显示 5更早的评论隐藏 5更早的评论

5 个评论
显示 3更早的评论隐藏 3更早的评论