シェルフフィルタ(ブーストフィルタ)

　バタワースフィルタだと減衰の傾きが-6dB/Octで落ちていくのだけど、ピンクノイズのような-3dB/Octの傾きはどうやって実現するのかというと、ハイシェルフフィルタ(ハイブーストフィルタ)を何段か通して上げてしまうらしい。周波数に比例して上げるのではなく棚壇上に上がるので、何段か通す必要がある。
　余りメジャーなフィルタではないのか、今まで持っていた参考書には載ってすらないし、検索で調べてもかかる量は少ない。それでも何とか出来たローシェルフフィルタは以下の感じ。グラフの上段が入力のインパルスで下段が周波数応答。赤線がインパルスで、青線がローシェルフフィルタの応答。

以下はscilabのソース。fcが棚壇上に上がる中央の周波数で、sが傾きのパラメータ。
fs=1000;
fc=10;
time=10;
gain=1;
s=100;
n=fs*time;
N=fs*time;

Implus=zeros(1,fs*time);
Implus(fs)=10;

w=(2*%pi*fc)/fs;
a=10^(gain/40);
fQ=sqrt(a+1/a)*(1/s-1)+2;
//fQ=1/sqrt(2);
p=sin(w)/(2*fQ);//?¿
c=2*sqrt(a)*p;
//c=sqrt(a/fQ)*sin(w);

B=;
B(1)=a*(a+1-(a-1)*cos(w)+c);
B(2)=2*a*(a-1-(a+1)*cos(w));
B(3)=a*(a+1-(a-1)*cos(w)-c);

A=;
A(1)=(a+1+(a-1)*cos(w)+c);
A(2)=-2*(a-1+(a+1)*cos(w));
A(3)=(a+1+(a-1)*cos(w)-c);

y=;
for i=1:N
if i==1
y(i)=B(1)*Implus(i)/A(1);
end

if i==2
y(i)=(B(1)*Implus(i)+B(2)*Implus(i-1)-A(2)*y(i-1))/A(1);
end

if 2;
FFTi=abs(fft(Implus',-1))*(2/fs);
FFTb=[];
FFTb=abs(fft(y,-1))*(2/fs);

Hz=(0:fs/n:(fs/n)*(n/2-1));
subplot(2,1,2)
plot2d(Hz(2:length(Hz)),FFTi(2:length(Hz)),logflag='ll',style=5)
plot2d(Hz(2:length(Hz)),FFTb(2:length(Hz)),logflag='ll',style=2)