Why does yulewalk produce an unexpected dip in the middle?

I think yulewalk is the right tool for the job here, but I can't figure out why it keeps producing an ugly dip in the middle of the response. I'm open to other solutions.

Goal: Invert/improve a transfer function response (loudspeaker measurement) to improve its impulse response in preparation for windowing.

% Get data
clc;clear;close;format compact; format short g;
Fs = 48000; Nyq = Fs/2;
url = 'https://subaligner.s3.us-east-2.amazonaws.com/TF/S1210_S1210+AC+WB+85-20k_-2.9dB.txt';
Ref = readtable(url,'ReadVariableNames',false); 
Ref(end,:) = []; % remove extra row
Ref.Properties.VariableNames = {'frequency','magnitude','phase','coherence'};
f = Ref.frequency;

% Convert to log spaced (this is the only way I've found to smooth logarithmically)
frequencyLog = round(logspace(log10(1),log10(Nyq),1000))';
frequencyLog = unique(frequencyLog);
rowsLog = ismember(f,frequencyLog);
TFlog = Ref(rowsLog,:);

% smooth
TFlogSmooth = smoothdata(TFlog.magnitude);

% interp back to linear
Ref.magnitudeSmooth = interp1(TFlog.frequency,TFlogSmooth,f,'pchip',0);

% find mean magnitude 50Hz to 18kHz
meanMag = mean(Ref.magnitudeSmooth(50:18000));

% mirror magnitude around mean
mean2zero = Ref.magnitudeSmooth - meanMag; % center around zero
Ref.magnitudeInv = mean2zero * -1; % mirror around zero

% limit magnitude to -20:+20
Ref.magnitudeInv(Ref.magnitudeInv < -20) = -20;
Ref.magnitudeInv(Ref.magnitudeInv > +20) = 20;

% design filter
fRad = f(1:Nyq) / Nyq;
fRad(1) = 0;
fRad(end) = 1;
m = rescale(db2mag(Ref.magnitudeInv(1:Nyq)));
order = 35; % ???
[b,a] = yulewalk(order,fRad,m);
% freqz(b,a,512)

% test filter
pulse = zeros(Fs,1);
pulse(Nyq) = 1;
yP = filter(b,a,pulse);
zP = fft(yP) / Fs;
magP = mag2db(abs(zP));
semilogx(f,magP + 110,f,Ref.magnitudeInv)
xlim([20 20000])

Answer 1

[Paul]

Hi Nathan,

Here's the original code from the question.

% Get data
clc;clear;close;format compact; format short g;
Fs = 48000; Nyq = Fs/2;
url = 'https://subaligner.s3.us-east-2.amazonaws.com/TF/S1210_S1210+AC+WB+85-20k_-2.9dB.txt';
Ref = readtable(url,'ReadVariableNames',false); 
Ref(end,:) = []; % remove extra row
Ref.Properties.VariableNames = {'frequency','magnitude','phase','coherence'};
f = Ref.frequency;
% Convert to log spaced (this is the only way I've found to smooth logarithmically)
frequencyLog = round(logspace(log10(1),log10(Nyq),1000))';
frequencyLog = unique(frequencyLog);
rowsLog = ismember(f,frequencyLog);
TFlog = Ref(rowsLog,:);
% smooth
TFlogSmooth = smoothdata(TFlog.magnitude);
% interp back to linear
Ref.magnitudeSmooth = interp1(TFlog.frequency,TFlogSmooth,f,'pchip',0);
% find mean magnitude 50Hz to 18kHz
meanMag = mean(Ref.magnitudeSmooth(50:18000));
% mirror magnitude around mean
mean2zero = Ref.magnitudeSmooth - meanMag; % center around zero
Ref.magnitudeInv = mean2zero * -1; % mirror around zero
% limit magnitude to -20:+20
Ref.magnitudeInv(Ref.magnitudeInv < -20) = -20;
Ref.magnitudeInv(Ref.magnitudeInv > +20) = 20;
% design filter
fRad = f(1:Nyq) / Nyq;
fRad(1) = 0;
fRad(end) = 1;
m = rescale(db2mag(Ref.magnitudeInv(1:Nyq)));

I noticed that m(1) was zero, so change it to 1.

m(1) = 1;
order = 35; % ???
[b,a] = yulewalk(order,fRad,m);

Frequency response of the filter.

[h,w] = freqz(b,a,512);
%figure;
%plot(w/pi,abs(h),fRad,m,f/Fs*2,rescale(db2mag(Ref.magnitudeInv))),grid
% test filter
pulse = zeros(Fs,1);

Changed the following line. Unit pulse response should start at n = 0.

%pulse(Nyq) = 1;
pulse(1) = 1;
yP = filter(b,a,pulse);
zP = fft(yP) / Fs;
magP = mag2db(abs(zP));
semilogx(f,magP + 110,f,Ref.magnitudeInv)
xlim([20 20000])

Make another plot comparing all of the responses: Fourier transform of yp, the design magnitude, the Ref.magnitudeInv (rescaled), and the frequency response of the filter.

figure
plot(f,abs(fft(yP)),fRad*Fs/2,m,f,rescale(db2mag(Ref.magnitudeInv)),w/pi*Fs/2,abs(h))
xlim([0 Fs/2])

We see that all the responses match pretty well for frequencies greater than 2000, but the frequency response of the filter doesn't really match the design magnitude at very low frequencies.

Zoom in

figure
plot(f,abs(fft(yP)),fRad*Fs/2,m,f,rescale(db2mag(Ref.magnitudeInv)),w/pi*Fs/2,abs(h),'o')
xlim([0 .2e4])

The filter response dips to nearly zero at 500, which then gets "magnified" when converted to dB.

My speculation is that the Yule-Walker filter will have trouble fitting that very narrow, nearly rectangular pulse for frequencies <200, which is a very, very small proportion of the overall frquency range.

[Nathan Lively]

@Paul Thank you so much for taking the time to make these thoughtful comments.

> I noticed that m(1) was zero, so change it to 1.

Thanks! I was totally backwards on my understanding of this. I thought it the IR you wanted to filter had to be centered.

> response of the filter doesn't really match the design magnitude at very low frequencies.

So what do you think are some solutions I could investigate?

I could use yulewalk on 200-20,000Hz. Then apply a low shelf to 0-200Hz.

Or is there another function more appropriate for this task than yulewalk?

I tried fir2, but soon realized that I didn't understand how to generate the minimum phase information from the magnitude response that I needed.

[Paul]

Sorry I can't be of more help w/o more information about what the problem actually is. I was actually making some educated guesses from the code.