Hi, How can I create an FIR from difference of two impulse responses?

Hello, I am trying to create an FIR and assocaited waveform for a trasnfer fucntion that is the magnitude difference between the two impulses responses. I have recorded both IR's but am struggling with getting the correction factor to be correct. Here is my code and results of the IRs and difference which is clearly wrong because it is reducing magntiude incorrectly. I can see kind of where it is going wrong which seems to be in taking the difference between them which results in the following magnitude graph and DBFS graph. I have inlcuded wav files as csv if trying to run.

when i mak apply 20*log(x) to each vector to find DBFS

% Define the sampling frequency
close all
%[x1,fs] = audioread('711_Tap1024.wav');
%[x2,fs] = audioread('5128_Tap1024.wav');
x1 = readmatrix('711_Tap1024.csv');
x2 = readmatrix('5128_Tap1024.csv');
fs = 48000;
Nfft = length(x1);
% Frequency vector (assume frequencies are normalized, or you need to normalize by fs/2)
f = linspace(0, fs, Nfft);
f = f';
% Define or import the frequency responses H1 and H2
H1 = fft(x1,Nfft);
H2 = fft(x2,Nfft);
% Compute the difference in frequency responses
H_diff = H1-H2;
% Inverse FFT to obtain the impulse response
y = real(ifft(H_diff,Nfft));
% Normalize the impulse response
y_new = 0.99*y/max(abs(y));
% Save the impulse response as a WAV file
audiowrite('impulse_response_diff.wav', y, fs);
% Plot the impulse response (optional)
figure;
plot(y)
y1=abs(fft(x1,Nfft));
y2=abs(fft(x2,Nfft));
y3=abs(fft(y,Nfft));
y1= 20*log10(y1);
y2= 20*log10(y2);
y3= 20*log10(y3);
figure
plot(f(1:Nfft/2), y1(1:Nfft/2));
hold on
plot(f(1:Nfft/2), y2(1:Nfft/2));
hold on
plot(f(1:Nfft/2), y3(1:Nfft/2));
xlim([20 20000])
xticks([100 1000 10000])
grid on
set(gca, 'XScale', 'log')
xlabel ('Frequency (Hz)');
ylabel ('Amplitude ');

[Pierce]

EDIT: I can get the corrections accuratly when i take the difference after abs(fft(x1)) and I am able to ahieve correct values for either power correciton or dB SPL correciton as vectors but by that point I am no longer able to ifft back to the time domain.

[Star Strider]

You need to subttract the complex results, not the abs results. Then, you should be able to invert the difference. (This is essentiially frequency-domain filtering, and while not recommended, is possible.)

That aside, it is difficult to understand what you want to do.

[Paul]

@Pierce

Can you show mathematically why the difference between H1 and H2 is needed?

From a comment below, it sounds like you have:

Y1 = H1*U

but you want to compute

Y2 = H2*U

based on knowing Y1, H2, and H1. If so, how would H2-H1 come into the analysis?

Answer 1

[Star Strider]

I slightly changed your original code here, in order to define my questions about it.

It is possible to use the firls function to approximatee the ‘y3’ plot. If you want to use it on your existing data, be miindful of the filter length (in this instance the filter order), since it should be less than one-third of the signal length.

This result is far from perfect, however it should provide a startiing point —

x1 = readmatrix('711_Tap1024.csv');
x2 = readmatrix('5128_Tap1024.csv');
fs = 48000;
Nfft = length(x1);
% Frequency vector (assume frequencies are normalized, or you need to normalize by fs/2)
f = linspace(0, fs, Nfft);
f = f';
% Define or import the frequency responses H1 and H2
H1 = fft(x1,Nfft);
H2 = fft(x2,Nfft);
% Compute the difference in frequency responses
H_diff = H1-H2;
% Inverse FFT to obtain the impulse response
y = H_diff;
% % y = real(ifft(H_diff,Nfft));
% Normalize the impulse response
y_new = 0.99*y/max(abs(y));
% Save the impulse response as a WAV file
% % audiowrite('impulse_response_diff.wav', y, fs);
% Plot the impulse response (optional)
% % figure;
% % plot(y)
y1=abs(fft(x1,Nfft));
y2=abs(fft(x2,Nfft));
y3=abs(fft(y,Nfft));
y1= 20*log10(y1);
y2= 20*log10(y2);
y3= 20*log10(y3);
figure
plot(f(1:Nfft/2), y1(1:Nfft/2), DisplayName='y_1');
hold on
plot(f(1:Nfft/2), y2(1:Nfft/2), DisplayName='y_2');
hold on
plot(f(1:Nfft/2), y3(1:Nfft/2), DisplayName='y_3');
xlim([20 20000])
xticks([100 1000 10000])
grid on
set(gca, 'XScale', 'log')
xlabel ('Frequency (Hz)');
ylabel ('Amplitude ');
legend('Location','best')

% ylim([-50 50])
format shortG
fn = fs/2

fn = 
       24000

Lv = f < fn;
fstop = (y3(Lv) < -275);
stopbands = f(fstop)

stopbands = 3×1

1.0e+00 * 4926.7 15343 20551

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

firn = 150;
fv = f < fn;
firf = f(fv);
a = y3(fv);
ftype = 'hilbert';
b = firls(firn, firf, a, ftype);

figure
freqz(b, 1, 2^16, fs)

.

[Pierce]

Hi Star Strider, thank you for responding!

Essentially what I am trying to do is create a correction curve for the gear that recorded H1. The two impulse responses are of a Hi-fidelity earbud recorded through a microphone that is either paird with an IEC Type 3.3 711 coupler or a B&K 5128 ear simulator. I am trying to get the difference between the two mic apparatices. Looking at my second graph the H2 is typicaly quieter below 1kHz and then is variable above it. I want to get the difference between these as a correction for H1 results so that whenever i take measurement on the 711 coupler, I apply my correction curve and see what the results would be like if i measure on the ear simulator. I can figure out the trasnfer function as a vector but not sure how to convert back to the time domain. I know this type of method is not alwasy ideal to use and just measuring on said gear is better, but it is a practice problem I am giving my self to get better aquainted with creating time domain FIRs.

One thought I had is apply the trasnfer funciton to a test sine sweep and then deconvolve with the orignal to ge the time domain response of the transfer function?

[Pierce]

So I am essentially trying to achieve this response in these graphs which are the same data just different units. dBFS vs power. but trying to get this as a result from taking the FFT of an FIR waveform.

[Star Strider]

My pleasure!

I am having problems understanding what you are doing. The filter based on ‘H_diff’ has a lowpass characteristic with a cutoff frequency of abour 18.5 kHz. It might be easier to just design a filter to do that, using the fir1 function. I did that as well here.

As I understand it, you want to do something liike this —

x1 = readmatrix('711_Tap1024.csv');
x2 = readmatrix('5128_Tap1024.csv');
fs = 48000;
Nfft = length(x1);
% Frequency vector (assume frequencies are normalized, or you need to normalize by fs/2)
f = linspace(0, fs, Nfft);
f = f';
% Define or import the frequency responses H1 and H2
H1 = fft(x1,Nfft);
H2 = fft(x2,Nfft);
% Compute the difference in frequency responses
H_diff = H1-H2;
fn = fs/2;
fp = f <= fn;                               % Logical Vector<

figure
plot(f(fp), mag2db(abs(H1(fp))*2), DisplayName='H_1');
hold on
plot(f(fp), mag2db(abs(H2(fp))*2), DisplayName='H_2');
hold on
plot(f(fp), mag2db(abs(H_diff(fp))*2), DisplayName='H_{diff}');
xlim([20 20000])
xticks([100 1000 10000])
grid on
set(gca, 'XScale', 'log')
xlabel ('Frequency (Hz)');
ylabel ('Amplitude ');
legend('Location','best')

format shortG
a = mag2db(abs(H_diff(fp))*2);
an = a - max(a);
idx = find(diff(sign(an + 6)));
hpf = f(idx)

hpf = 6×1

1.0e+00 * 5818.2 12434 14311 17736 17783 18534

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

figure
plot(f(fp), mag2db(abs(H_diff(fp))*2))
hold on
plot(f(fp), an)
hold off
grid
yline(-6, '--k', '-6 dB')
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)')
title('H_{diff}')

fn = fs/2;

firn = fix(numel(x1)/4);
firf = f(fp)/fn;
% a = mag2db(abs(H_diff(fp)));
% a = a - max(a);
ftype = 'hilbert';
b1 = firls(firn, firf, an, ftype);

figure
freqz(b1, 1, 2^16, fs)

b1 = firls(firn, firf, an);

figure
freqz(b1, 1, 2^16, fs)

b2 = fir1(firn, max(hpf)/fn);

figure
freqz(b2, 1, 2^16, fs)

% set(subplot(2,1,1), 'XScale','log', 'XLim',[20 20000])
% set(subplot(2,1,2), 'XScale','log', 'XLim',[20 20000])

.

[Pierce]

Hi Star,

Thank you again for answering. this is not quite what I am looking for and have updates since. Essentially I have 2 impulses responses that give me a SPL frequency response in dB after abs fft and converting to decibel scale. I am trying to find the deltas between them as an FIR so that when i apply the FIR to IR 1 it shows me the FR of IR2. so lets say at 100Hz IR2 is 5dB quiter than IR1. that means my resulting FIR should have a -5dB value at 100Hz so thaqt in the future when I apply the FIR to IR1 I will get the result as if I meaured form the IR2 system. hopefully that makes sense.

I have since found I need to divide H2./H1 to get the right data. I am now getting the right results, HOWEVER my time domain plot of the correciton curve does not look correct to me even though when I take the FFT of it I get the result i want.

Here is updated code.

% Define the sampling frequency

close all

%[x1,fs] = audioread('711_Tap1024.wav');

%x2,fs] = audioread('5128_Tap1024.wav');

x1 = readmatrix('711_Tap1024.csv');

x2 = readmatrix('5128_Tap1024.csv');

fs = 48000;

Nfft = length(x1);

% Frequency vector (assume frequencies are normalized, or you need to normalize by fs/2)

f = linspace(0, fs, Nfft);

f = f';

%define Target Curve

x1FR = 20*log10(abs(fft(x1,Nfft)));

x2FR = 20*log10(abs(fft(x2,Nfft)));

targetCurve = x1FR-x2FR;

targetCurve = targetCurve*-1;

% Define or import the frequency responses H1 and H2

% H1 = abs(fft(x1,Nfft));

% H2 = abs(fft(x2,Nfft));

H1 = fft(x1,Nfft);

H2 = fft(x2,Nfft);

% Compute the difference in frequency responses divide

%H_diff = H1-H2;

H_diff = H2./H1;

% Inverse FFT to obtain the impulse response

y = real(ifft(H_diff,Nfft));

% Normalize the impulse response

y = 0.99*y/max(abs(y));

% Save the impulse response as a WAV file

audiowrite('impulse_response_diff.wav', y, fs);

% Plot the impulse response (optional)

figure;

plot(y)

% plot result of y

y_new=abs(fft(y,Nfft));

y_new= 20*log10(y_new);

figure

plot(f(1:Nfft/2), x1FR(1:Nfft/2));

hold on

plot(f(1:Nfft/2), x2FR(1:Nfft/2));

hold on

plot(f(1:Nfft/2), targetCurve(1:Nfft/2));

hold on

plot(f(1:Nfft/2), y_new(1:Nfft/2));

xlim([20 20000])

xticks([100 1000 10000])

grid on

set(gca, 'XScale', 'log')

xlabel ('Frequency (Hz)');

ylabel ('Magnitude');

legend('H1', 'H2', 'Target Curve','Correction Curve')

[Star Strider]

I cannot get a good approximation to your desired transfer function using the Signal Proceessing Toolbox functions.

Note that there was some high-freequency noise at the high end of the ‘H_diff’ transfer functiion. I simply eliminated that by truncating the vectors.

Using the System Identification Toolbox functions, I was able to get a good approximation in the continuous s domain. I could not determine a sampling frequency from your available data, so you will have to transform the estimated transfer function to a discrete (digital) filter using the c2d function. It should work, however I can offer no guarantees. When you get the discrete numerator and denominator vectors, it would be best to use the tf2sos function, and then use the second-order-section representatiion to filter your signal, using the filtfilt function.

% Define the sampling frequency
% close all
%[x1,fs] = audioread('711_Tap1024.wav');
%x2,fs] = audioread('5128_Tap1024.wav');
x1 = readmatrix('711_Tap1024.csv');
x2 = readmatrix('5128_Tap1024.csv');
fs = 48000;
Nfft = length(x1);
% Frequency vector (assume frequencies are normalized, or you need to normalize by fs/2)
f = linspace(0, fs, Nfft);
f = f';
%define Target Curve
x1FR = 20*log10(abs(fft(x1,Nfft)));
x2FR = 20*log10(abs(fft(x2,Nfft)));
targetCurve = x1FR-x2FR;
targetCurve = targetCurve*-1;
% Define or import the frequency responses H1 and H2
% H1 = abs(fft(x1,Nfft));
% H2 = abs(fft(x2,Nfft));
H1 = fft(x1,Nfft);
H2 = fft(x2,Nfft);
% Compute the difference in frequency responses divide
%H_diff = H1-H2;
H_diff = H2./H1;
% Inverse FFT to obtain the impulse response
y = real(ifft(H_diff,Nfft));
% Normalize the impulse response
y = 0.99*y/max(abs(y));
% Save the impulse response as a WAV file
audiowrite('impulse_response_diff.wav', y, fs);
% Plot the impulse response (optional)
figure;
plot(y)

% plot result of y
y_new=abs(fft(y,Nfft));
y_new= 20*log10(y_new);
figure
plot(f(1:Nfft/2), x1FR(1:Nfft/2));
hold on
plot(f(1:Nfft/2), x2FR(1:Nfft/2));
hold on
plot(f(1:Nfft/2), targetCurve(1:Nfft/2), LineWidth=3);
hold on
plot(f(1:Nfft/2), y_new(1:Nfft/2));
xlim([20 20000])
xticks([100 1000 10000])
grid on
set(gca, 'XScale', 'log')
xlabel ('Frequency (Hz)');
ylabel ('Magnitude');
legend('H1', 'H2', 'Target Curve','Correction Curve', Location='best')

fn = fs/2;
fp = f(1:Nfft/2);
H_diff = H_diff(1:Nfft/2);
Lv = fp<=2.13E+4;
fp = fp(Lv);
H_diff = H_diff(Lv);

nnz(Lv)

ans = 454

wp = fp/fn * pi;

% H_diff = H_diff(Lv);
% H_diff(1:10)
H_diff(1) = H_diff(2);
H_diffn = abs(H_diff) / max(abs(H_diff));

[b,a] = invfreqz(wp, abs(H_diffn(1:numel(fp))), 64, 64);

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.705863e-21.

[hfz1,wfz1] = freqz(b, a, 2^16);
b = firls(64, wp, abs(H_diffn(1:numel(wp))), 'hilbert')

b = 1×65

-0.0005 -0.0007 -0.0012 -0.0007 -0.0010 -0.0009 -0.0008 -0.0009 -0.0007 -0.0012 -0.0010 -0.0012 -0.0018 -0.0013 -0.0020 -0.0018 -0.0017 -0.0021 -0.0012 -0.0013 -0.0023 -0.0020 -0.0008 -0.0029 -0.0000 -0.0147 0.0116 -0.0168 0.0138 -0.0799

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[hfz2,wfz2] = freqz(b, 1, 2^16);

frd = idfrd(H_diffn(1:numel(fp)), fp, 0)

frd =
IDFRD model.
Contains Frequency Response Data for 1 output(s) and 1 input(s).
Response data is available at 454 frequency points, ranging from 0 rad/s to 2.126e+04 rad/s.
 
Status:                                                         
Created by direct construction or transformation. Not estimated.
 

sys = tfest(frd, 6)

Warning: Phase information was added to the real and non-negative signals. See "addMinPhase" command for more information.

Warning: Estimated model has a lower order than requested.

sys =
 
  0.1372 s^5 + 2690 s^4 + 4.837e07 s^3 + 2.67e11 s^2 + 2.538e15 s + 5.58e17
  -------------------------------------------------------------------------
    s^5 + 3918 s^4 + 1.918e08 s^3 + 3.919e11 s^2 + 7.102e15 s + 2.604e18
 
Continuous-time identified transfer function.

Parameterization:
   Number of poles: 5   Number of zeros: 5
   Number of free coefficients: 12
   Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                                
Estimated using TFEST on frequency response data "frd".
Fit to estimation data: 93.77%                         
FPE: 0.0004776, MSE: 0.000455                          
 

format longE
b = sys.Numerator

b = 1×6

1.0e+00 * 1.371929903001000e-01 2.690026894243068e+03 4.837384115249185e+07 2.669725032767512e+11 2.537721438997358e+15 5.580319150849004e+17

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

a = sys.Denominator

a = 1×6

1.0e+00 * 1.000000000000000e+00 3.918119344768781e+03 1.917799060958707e+08 3.918664071344840e+11 7.102028871891668e+15 2.603730705315592e+18

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

format short

figure
compare(frd, sys)

figure
tiledlayout(3,1)
nexttile
plot(wp, abs(H_diffn(1:numel(fp))))
xlim([0 pi])
title('Original Transfer Function')
% hold on
% plot(fp, ftf)
% xline(2.14E+4,'--k')
% hold off
grid
nexttile
plot(wfz1, abs(hfz1))
title('‘invfreqz’ Estimate')
grid
nexttile
plot(wfz2, abs(hfz2))
title('‘firls’ Estimate')
grid

This is the best that I can do with your data. This is definitely not straightforward.

.