Frequency Response to IR calculation yields shallower Magnitude level?

5 次查看（过去 30 天）
Michaela Warnecke
Michaela Warnecke 2023-10-4
评论： Mathieu NOE 2023-10-9
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?
  2 个评论
Paul
Paul 2023-10-6
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
Michaela Warnecke 2023-10-9
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:

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Michaela Warnecke
Michaela Warnecke 2023-10-9
See above response to Paul's note.

更多回答（1 个）

Mathieu NOE
Mathieu NOE 2023-10-6
Ran in:
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 (°)');

类别

在 Help Center 和 File Exchange 中查找有关 Filter Design 的更多信息

标签

产品


版本

R2022b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by