Why does yulewalk produce an unexpected dip in the middle?
1 次查看(过去 30 天)
显示 更早的评论
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1093890/image.png)
% 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])
0 个评论
采纳的回答
Paul
2022-8-11
编辑:Paul
2022-8-11
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 个评论
Paul
2022-8-11
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.
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Multirate Signal Processing 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!