Frequency Response to IR calculation yields shallower Magnitude level?

Hi all,

I have a frequency response curve of an audio recording. All I know is that is spans some range (see below) and was (most likely) sampled at fs = 48k.

I need to convert it to an IR to apply some specific filtering and am using the following code for that conversion:

FFT_length = 1024;
fs = 48e3;
%freq_axis = 1xN vector containing the range of frequencies in Hz ([71:20e3] in some octave steps)
%FR_data = 1xN vectora containing the magnitde values at each freq in freq_axis

%convert frequencies to normalized values to generate IR filters
filter_freqs = [0 freq_axis./(fs/2) 1];
%Create time-domain IRs
IR_sig = fir2(FFT_length, filter_freqs, [FR_data(1)-1 FR_data FR_data(end)+1]);

This creates an IR signal, which I then convert back into the frequency domain to validate against the original FR using the following:

cutoff_freq = fs/2;
HRTF_mag = abs(fft(IR_sig, FFT_length));

%Return only  from 0 to pi
Mag = 20*log10(HRTF_mag(1:(FFT_length/2)+1));
f_vec = [0:(cutoff_freq*2)/FFT_length:cutoff_freq];

When I plot the original Frequency response curve (FR_data, left) and also plot the IR-converted, frequency domain signal (Mag, right), I can see a close match in the frequency domain, but the magnitude is entirely off. You can see that the original data (left) shows a difference of ~30 dB between the notch at 2.5k and the peak at 4.1k. The same notch and peak are only about 5dB apart in the deduced FR (Mag).

What am I missing/doing wrong?

[Paul]

Hi Michaela,

You're more likely to get help if you upload the data needed to run your code and recreate your results. Save freq_axis and FR_data (if that's all that's needed) to a .mat file and add the file to your question using the paperclip icon in the Insert ribbon. Also, can you add the plotting commands used to make the plots.

[Michaela Warnecke]

Hi Paul,

Thanks for your note. I have since figured out that the issue was that values going into fir2 are dB and need to be linear, so the fix is to do the following for this case:

b = fir2(n, f, 10.^(m./20)) 

About adding data: It's possible to just run this with made-up data and show the issue:

FFT_length = 1024;
fs = 48e3;

%making up FR_data and freq_axis
FR_data =  30 + (100-30)*rand(1,100);% selecting min/max based on figure above. FR_data
% containing the magnitde values at each freq in freq_axis
freq_axis = linspace(70, 16e3, numel(freq_axis)); % selecting min/max based on figure above. freq_axis
% containing the range of frequencies in Hz ([71:20e3] in some octave steps)

%convert frequencies to normalized values to generate IR filters
filter_freqs = [0 freq_axis./(fs/2) 1];
%Create time-domain IRs
IR_sig = fir2(FFT_length, filter_freqs, [FR_data(1)-1 FR_data FR_data(end)+1]);
cutoff_freq = fs/2;
HRTF_mag = abs(fft(IR_sig, FFT_length));
%Return only  from 0 to pi
Mag = 20*log10(HRTF_mag(1:(FFT_length/2)+1));
f_vec = [0:(cutoff_freq*2)/FFT_length:cutoff_freq];

%plot first data to show the same issue as in original ask
figure; set(gcf, 'pos', [80 400 1250 480])
subplot(121); semilogx(freq_axis, FR_data, 'k'); xlim([40 21e3]); grid on
subplot(122); semilogx(f_vec, Mag, 'b'); xlim([40 21e3]); grid on

%fix issue by ajusting for dB/linear
IR_sig2 = fir2(FFT_length, filter_freqs, 10.^([FR_data(1)-1 FR_data FR_data(end)+1]./20));
HRTF_mag2 = abs(fft(IR_sig2, FFT_length));
%Return only  from 0 to pi
Mag2 = 20*log10(HRTF_mag2(1:(FFT_length/2)+1));
f_vec2 = [0:(cutoff_freq*2)/FFT_length:cutoff_freq];

%plot updated data
subplot(121); hold on; semilogx(f_vec2, Mag2, 'r'); xlim([40 21e3])
text(80, 90, 'FR_data', 'Color', 'k', 'Interpreter','none', 'FontSize', 14)
text(80, 88, 'Mag2', 'Color', 'r', 'Interpreter','none', 'FontSize', 14)
subplot(122); hold on; semilogx(f_vec2, Mag2, 'r'); xlim([40 21e3])
text(80, 90, 'Mag', 'Color', 'b', 'Interpreter','none', 'FontSize', 14)
text(80, 88, 'Mag2', 'Color', 'r', 'Interpreter','none', 'FontSize', 14)

This produces the following:

Answer 1

[Michaela Warnecke]

See above response to Paul's note.

Answer 2

[Mathieu NOE]

hi

I have done something similar in the past - see below FYI

I start with a given FIR filter then computation of the transfer function and backward to the FIR ( = impulse response of a TF)

% 
Fs = 1e3;
coef_fs = 2; % 2 for Nyquist (exact match) , measurement factor can be 2.56 
Freq = linspace(0,Fs/coef_fs,200); % make sure freq vector goes up to Fs/2

b = fir1(48,[0.2 0.3]);   % Window-based FIR filter design
frf = freqz(b,1,Freq,Fs);

%% FIR obtained with ifft method

if mod(length(frf),2)==0 % iseven
    frf_sym = conj(frf(end:-1:2));
else
    frf_sym = conj(frf(end-1:-1:2));
end
    
fir = real(ifft([frf frf_sym]));
% resample fir (impulse response) if frf data is not provided up to Nyquits frequency (Fs/2)
corfactor = 2/coef_fs;
x_og = (1:length(fir));
x_cor = linspace(1,length(fir)*corfactor,length(fir));
fir = interp1(x_og,fir,x_cor);
fir = fir(1:50); % truncation is possible if fir decays enough
frfid = freqz(fir,1,Freq,Fs);

% alternative with invfreqz
ITER = 1e3;
% FIR filter design
NB = 50; % 
NA = 0;
W = 2*pi*Freq/Fs;
Wt = ones(size(W));
TOL = 1e-2;
[fir2,A] = invfreqz(frf,W,NB,NA,Wt,ITER,TOL);
frfid2 = freqz(fir2,1,Freq,Fs);


figure(1)
subplot(311),plot(1:length(b),b,'-*',1:length(fir),fir,1:length(fir2),fir2);
legend('FIR model','identified FIR (fft)','identified FIR (invfreqz)');
xlabel('samples');
ylabel('amplitude');
subplot(312),plot(Freq,20*log10(abs(frf)),Freq,20*log10(abs(frfid)),Freq,20*log10(abs(frfid2)),'*');
legend('input model','identified FIR model (fft)','identified FIR model (invfreqz)');
xlabel('frequency (Hz)');
ylabel('FRF modulus (dB)');
subplot(313),plot(Freq,180/pi*angle(frf),Freq,180/pi*angle(frfid),Freq,180/pi*angle(frfid),'*');
legend('input model','identified FIR model (fft)','identified FIR model (invfreqz)');
xlabel('frequency (Hz)');
ylabel('FRF angle (°)');

[Michaela Warnecke]

Thanks, Mathieu. Nice solution. I figured out that my issue was a simple dB/linear misrepresentation (see above).

[Mathieu NOE]

ok , good to hear

all the best for the future