2
$\begingroup$

I am trying to find the reverberation time of a room using the Schroeder method (i.e., Reverse-time integration method). Therefore, impulse responses should be measured first.

There are many ways to obtain the impulse response (IR) of a room, such as maximum length sequence, exponential sine sweep method, etc. First, I tried to obtain IR from frequency response performing IFFT, because I can easily obtain the frequency response using measuring devices and software (LMS Test.Lab). The schematic diagram of the measuring equipment is as follows:

enter image description here

Here, Input is the signal from the Amp. and output is the sound pressure measured at the Microphone. Then I get the frequency response as follows: enter image description here

From that, I try to obtain the filtered impulse response (1/3octave band) using MATLAB code. Below is an excerpt from some of the matlab code to implement this.

%%  Design 1/3 octave Band pass filter (Butterworth )
% yf: single-sided spectrum from LMS Test.Lab
% Nf: number of frequency response data
% Nt: number of time data, usually Nt=2*(Nf-1)

%% Setting parameters
fcMin = 100; % Lower center freq.
fcMax = 8000; % Upper center freq. 
bw = 1/3; 
octs = log2(fcMax/fcMin);
bmax = ceil(octs/bw);
fc = fcMin*2.^( (0:bmax) * bw ); % center frequencies
fl = fc*2^(-bw/2); % lower cutoffs
fu = fc*2^(+bw/2); % upper cutoffs
numBands = length(fc);
b = cell(numBands,1);
a = cell(numBands,1);


frfmat=zeros(numBands,Nf); 
htmat=zeros(numBands,Nt); 
 
%% Butterworth Band pass filter

 for nn = 1:length(fc)

    [b{nn},a{nn}] = butter(3, [fl(nn) fu(nn)]/(fs/2), 'bandpass');
    [h,f]=freqz(b{nn},a{nn},Nf,fs);
    
    frfmat(nn,:)=(yf.') .* abs(h.');  % Multiplying amplitude of filter response and single-sided spectrum in narrow band.
    htmat(nn,:) = (IFFT_hy(frfmat(nn,:),Nf)).'; % synthesize impulse response from single-sided spectrum;
    
     
end

Obviously, the frequency response (yf) is multiplied by the amplitude of the Butterworth filter (h) to preserve the phase information of the frequency response, and then inverse FFT (IFFT_hy) is used to synthesize the filtered impulse response.

Is this the correct way of obtaining the impulse response filtered by 1/3 octave band, or is anything that needs to be corrected?

------------------------ Additional question ----------------------------

Based on Hirmar's constructive feedback, 1/3 octave band filtering was performed again using butterworth filter.

Three cases were used. The first method is to perform inverse FFT after multiplying the frequency response by the amplitude of the butterworth filter in the frequency domain as described above.

The second and third methods are to apply zero-phase Butterworth filter(filtfilt function) and Butterworth filter(filter function) to time data (yt), respectively.

Below is the pseudo MATLAB code to obtain the above-mentioned results.

for nn = 1:length(fc)

    [b{nn},a{nn}] = butter(3, [fl(nn) fu(nn)]/(fs/2), 'bandpass');
    
    % Applying filter in Frequency domain 
    [h,f]=freqz(b{nn},a{nn},Nf,fs);
    frfmat(nn,:)=(yf.') .* abs(h.');  % First case!!
    htmat_fd(nn,:) = (IFFT_DSS_hy(frfmat(nn,:),Nf)).'; % synthesize impulse response from single-sided spectrum;
    
   
    % Applying filter in time domain 
    htmat_filtfilt(nn,:) = filtfilt(b{nn},a{nn},yt) ; % Second case!!
    htmat_filter(nn,:) = filter(b{nn},a{nn},yt) ; % Third case!!
    

end

I wanted to upload the calculated result, but the function to upload an image using Imgur is not currently available. It would be appreciated if you could tell us the part that needs to be improved in the process of filtering with a 1/3 Octave band using the Butterworth filter based on the second method.

$\endgroup$
2
  • 1
    $\begingroup$ If you need to preserve the phase then a Butterworth (IIR) doesn't seem like a logical choice but, I don't know the method. $\endgroup$ Nov 2, 2022 at 8:46
  • $\begingroup$ @aconcernedcitizen You'r right, the Bessel filter is the best Linear Phase filter in the analog domain. I have used them. Really neat to see on a scope: en.wikipedia.org/wiki/Bessel_filter . But since I never needed a flat gain response, an adjustment might need to be made. $\endgroup$
    – rrogers
    Nov 8, 2022 at 20:01

1 Answer 1

0
$\begingroup$

Is this the correct way of obtaining the impulse response filtered by 1/3 octave band,

No. This will not work at all. Frequency domain processing is mathematically quite complicated and I would not recommend it unless you have a thorough understanding of the underlying mathematical fundamentals (circular convolution, time-domain aliasing, windowing, etc.)

How to obtain filtered impulse response from frequency response?

In your case you can calculate broad band impulse response first and then apply the third octave filters in the time domain.

Your data needs to be on a complete FFT frequency grid and contain both magnitude and phase with decent signal to noise ratio at all frequencies. IF that's the case, you can simply make the spectrum conjugate symmetric and apply an inverse FFT. Than inspect the resulting impulse response for obvious errors: pre-ringing, pre-echoes, excessive noise such as very low frequency noise, DC drift, line hum (50/60Hz and harmonics), etc. If it looks good just apply the BW filters one at a time.

If your spectrum data isn't in the right format, than things get a lot more complicated.

$\endgroup$
4
  • $\begingroup$ Dear Hilmar, thank you again for your detailed instruction. I've revised the text based on your feedback. I would really appreciate it if you could review it when you have time! $\endgroup$
    – Toney
    Nov 3, 2022 at 1:48
  • $\begingroup$ As you said, the time data obtained using the Inverse FFT shows 50-60 Hz harmonics noise components possibly due to the line hum. If I use a bandpass filter to find the 1/3 octave band impulse response above 100 Hz, will these components still be a problem? $\endgroup$
    – Toney
    Nov 3, 2022 at 2:27
  • $\begingroup$ filter() is definitely the way to go. All other methods are non-causal which you really do NOT want here. Whether the noise is a problem or not depends on the details. I can't tell without looking at the actual data. You can find specification for 3rd octave filters in IEC 61260-1:2014 or ANSI S1-6-2016 . I don't know out of the top of my hat whether 3rd BW bandpass generally qualifies or not. That depends also on the "class" of 3rd octave filter you need. $\endgroup$
    – Hilmar
    Nov 3, 2022 at 13:46
  • $\begingroup$ I am grateful to you for sharing your experience. I will study more on this topic. Have a happy weekend! $\endgroup$
    – Toney
    Nov 4, 2022 at 1:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.