Я тож посидул покумекал, Гилберт там тож через ффт делается, так что результат такой же
Только вот импульс перед прогонкой через эту лабуду надо увеличивать раз в 8, иначе точность падает - критично на пример, на анти-алиасигн FIR-фильтре для ресемплера.
Вот окончательный вариант
function y=mpestimate(b,N,s)
% b - impulse
% N - impulse points, use ~ 2^nextpow2(length(b))*8
% s - output impulse length
mag=abs(fft(b,N));
% clip to avoid zeros ( log(0)=-Inf )
cutoff=-100; %db
thr=max(mag).*10^(cutoff/20);
toosmall=mag<thr;
mag(toosmall)=thr;
sig(1:(N/2))=sign(linspace(1,(N/2),(N/2)));
sig((N/2)+1)=0;
sig((N/2)+2:N)=sign(linspace(-1,-(N/2)-1,(N/2)-1));
sig(1)=0;
ceps=real(ifft(log(mag)));
ph=-1i*fft(sig'.*ceps);
rec=mag.*exp(1i*ph);
recu=ifft(rec);
y=recu(1:s);
end
По поводу клиппера - оно то конечно так, только на практике не важно ибо низкие частоты на входе сами смещают рабочую точку ограничителя
Мы ж не чистый синус туда подаем, на выходе куча четных и нечетных гармоник + интермодуляция, вообще входным фильтром все рулится...
Единственное, что можно добавить, так эту имитацию выпрямителя, а-ля катодный повторитель в меса-ректо
float _a01,_b00, _a11,_b10,_yg1,_yk1;
float y1;
inline float process(float x){
float y;
y=y1*_yg1-_yk1;
if(x>y){
y=x;
y1=-y1*_a01 +x*_b00;
}else{
y1=-y1*_a11 -1*_b10;
y=y1*_yg1-_yk1;
}
return y;
}
Коэфициенты вычисляем, как для обычного IIR первого порядка
void setFreq(float foff,float fon){ // f=f/fs
_a01=2*fon*pi-1;
_a11=2*foff*pi-1;
_b00=2*fon*pi;
_b10=2*foff*pi;
_yg1=(-foff+fon)/fon;
_yk1=foff/fon;
}Где foff,fon - частота(временная постоянная) релиза/атаки выпрямителя, в герцах/Fs
ps. Вникать в теорию и писать на C++ нужно, если хочется в живую полабать, а не гонять готовое через матлаб... хотя к приведению к МФ это не касается
pps. Спектр после деконволюци(в том числе и через pwelch) еще хорошо в логарифмическую шкалу по оси частот перевести, потом через медианный фильтр прогонять, а потом обратно в линейную... красивая картинка получается (на слух тоже). Это если референс не спец сигнал(свип, на пример), а кусок каккой-то гитарной партии..
Размещение рекламы