# Why does yulewalk produce an unexpected dip in the middle?

1 vue (au cours des 30 derniers jours)
Nathan Lively le 10 Août 2022
Commenté : Paul le 11 Août 2022
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(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
m = rescale(db2mag(Ref.magnitudeInv(1:Nyq)));
order = 35; % ???
% 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])
##### 0 commentairesAfficher -2 commentaires plus anciensMasquer -2 commentaires plus anciens

Connectez-vous pour commenter.

### Réponse acceptée

Paul le 11 Août 2022
Modifié(e) : Paul le 11 Août 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(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
m = rescale(db2mag(Ref.magnitudeInv(1:Nyq)));
I noticed that m(1) was zero, so change it to 1.
m(1) = 1;
order = 35; % ???
Frequency response of the filter.
[h,w] = freqz(b,a,512);
%figure;
% 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
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
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 commentairesAfficher AucuneMasquer Aucune
Nathan Lively le 11 Août 2022
Modifié(e) : Nathan Lively le 11 Août 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 le 11 Août 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.

Connectez-vous pour commenter.

### Catégories

En savoir plus sur Smoothing and Denoising dans Help Center et File Exchange

R2021a

### Community Treasure Hunt

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

Start Hunting!

Translated by