Why does yulewalk produce an unexpected dip in the middle?

1 view (last 30 days)
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])

Accepted Answer

Paul
Paul on 11 Aug 2022
Edited: Paul on 11 Aug 2022
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.
  2 Comments
Nathan Lively
Nathan Lively on 11 Aug 2022
Edited: Nathan Lively on 11 Aug 2022
@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
Paul on 11 Aug 2022
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.

Sign in to comment.

More Answers (0)

Products


Release

R2021a

Community Treasure Hunt

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

Start Hunting!