Properties of smoothing filters
samplingF = 1024;
stepT = 1/samplingF;
seqL = 2;
filterL = 5;
start = 100;
controlVar(start,100, filterL,10);
contF(x) = if x < start * stepT then 0 else
if x < ( start + filterL ) * stepT then 1/integer(filterL) else 0;
plot(contF,
(start - 20) * stepT, 2,
(start + filterL + 20) * stepT, -2 );
f(x)=sampleL(contF,0,seqL,seqL*samplingF);
plot(f);
ff(x)=fft1(f);
ffn(x)=ff(x)/maxS(ff); //normalize, or multiply by the samplingF?
ffp(x)=abs(ffn(x));
plot(ffp); //autoScale, should I divide by maxS(ffp)?
db(x)=8.685890*log(ffp(x));
plot(db);
iff(x)=ifft1(ff);
plot(iff);
//new stuff
flr(x)=firResp(fl); //what is fl?
plot(flr,-.1,1.1,Pi+.1,-.5);
Compare different window functions
samplingF = 2048; //Hz
samplingTime = 1; //seconds
filterLength = 41; //points
shape(x) = if x < ( samplingF / 2 ) / 3 then 1 else 0; //say limit to a third of the interval
fr(x)=sampleN( shape, 0, 1, samplingF/2+1 ); //freq. response
//plot(fr);
fk(x)=ifft1(fr); //filter kernel
//plot(fk);
fktr(x)=truncateS(shrotS(fk,filterLength/2),0,filterLength-1);
fkn(x)=fktr(x)/sumS(fktr); //normalized kernel
//plot(fkn);
hammingW(x) = 0.54 - 0.46 * cos( 2 * Pi * x / filterLength ); //Hamming
blackmanW(x) = 0.42 - 0.5 * cos( 2 * Pi * x / filterLength ) + 0.08 * cos( 4 * Pi * x / filterLength ); //Blackman
wfkh(x)=hammingW(x/stepS(fkn))*fkn(x); //the Hamming windowed filter
wfkb(x)=blackmanW(x/stepS(fkn))*fkn(x); //the Blackman windowed filter
//check the resulting filters frequency response
efr(x)=firResp(fkn); //estimated freq. response non-windowed
efrh(x)=firResp(wfkh); //estimated freq. response Hamming
efrb(x)=firResp(wfkb); //estimated freq. response Blackman
plot(efr,efrh,efrb,-.1,1.1,Pi+.1,-.5);
db(x)=8.685890*log(x);
efrdb(x)=db(efr(x)); //same as above but in db
efrdbH(x)=db(efrh(x));
efrdbB(x)=db(efrb(x));
plot(efrdb,efrdbH,efrdbB,-.1,1,Pi+.1,-70);
//alternatively do the same through the fourier transform of the filter kernel
fknp(x)=padS(fkn,2^nextpow2(samplingF*samplingTime)); //pad to the proper length
fkn1f(x)=fft1(fknp); //take fft
fkn1a(x)=abs(fkn1f(x))*samplingF/2; //take abs and normalize
fkn1db(x)=db(fkn1a(x)); //calculate in db
wfkp1(x)=padS(wfkh,2^nextpow2(samplingF*samplingTime)); //same for the Hamming windowed kernel...
wfkh1f(x)=fft1(wfkp1);
wfkh1a(x)=abs(wfkh1f(x))*samplingF/2;
wfkh1db(x)=db(wfkh1a(x));
wfkbp(x)=padS(wfkb,2^nextpow2(samplingF*samplingTime)); //same for the Blackman windowed kernel...
wfkb1f(x)=fft1(wfkbp);
wfkb1a(x)=abs(wfkb1f(x))*samplingF/2;
wfkb1db(x)=db(wfkb1a(x));
//plot(fkn1,wfkh1,wfkb1);
plot(fkn1a,wfkh1a,wfkb1a);
plot(fkn1db,wfkh1db,wfkb1db);
Compare filters with different lengths
samplingF = 2048; //Hz
samplingTime = 1; //seconds
filterLength1 = 41; //points
filterLength2 = 81; //points
shape(x) = if x < ( samplingF / 2 ) / 3 then 1 else 0; //say limit to a third of the interval
fr(x)=sampleN( shape, 0, 1, samplingF/2+1 ); //freq. response
//plot(fr);
fk1(x)=ifft1(fr); //filter kernel 1
fk2(x)=cloneS(fk1); //filter kernel 2
//plot(fk);
fkshr1(x)=shrotS(fk1,filterLength1/2);
fktr1(x)=truncateS(fkshr1,0,filterLength1-1);
fkn1(x)=fktr1(x)/sumS(fktr1); //normalized kernel 1
fkshr2(x)=shrotS(fk2,filterLength2/2);
fktr2(x)=truncateS(fkshr2,0,filterLength2-1);
fkn2(x)=fktr2(x)/sumS(fktr2); //normalized kernel 2
blackmanW(x,filterLength) = 0.42 - 0.5 * cos( 2 * Pi * x / filterLength ) + 0.08 * cos( 4 * Pi * x / filterLength ); //Blackman
wfk1(x)=blackmanW(x/stepS(fkn1),filterLength1)*fkn1(x); //the Blackman windowed filter 1
wfk2(x)=blackmanW(x/stepS(fkn2),filterLength2)*fkn2(x); //the Blackman windowed filter 2
plot(wfk1,wfk2);
//check the resulting filters frequency response
efr1(x)=firResp(wfk1); //estimated freq. response length 1
efr2(x)=firResp(wfk2); //estimated freq. response length 2
plot(efr1,efr2,-.1,1.1,Pi+.1,-.5);
db(x)=8.685890*log(x);
efrdb1(x)=db(efr1(x)); //same as above but in db 1
efrdb2(x)=db(efr2(x)); //same as above but in db 2
plot(efrdb1,efrdb2,-.1,1,Pi+.1,-70);
Filter inverse, fiter reverse
samplingF = 2048; //Hz
samplingTime = 1; //seconds
filterLength = 81; //points
shape(x) = if x < ( samplingF / 2 ) / 3 then 1 else 0; //say limit to a third of the interval
fr(x)=sampleN( shape, 0, 1, samplingF/2+1 ); //freq. response
//plot(fr);
fk(x)=ifft1(fr); //filter kernel 1
//plot(fk);
fkshr(x)=shrotS(fk,filterLength/2);
fktr(x)=truncateS(fkshr,0,filterLength-1);
fkn(x)=fktr(x)/sumS(fktr); //normalized kernel
blackmanW(x,filterLength) = 0.42 - 0.5 * cos( 2 * Pi * x / filterLength ) + 0.08 * cos( 4 * Pi * x / filterLength ); //Blackman
wfk(x)=blackmanW(x/stepS(fkn),filterLength)*fkn(x); //the Blackman windowed filter
fkInv(x)=inverseFilter(wfk); //inverse
//fkInv(x)=reverseFilter(wfk); //reverse
plot(wfk,fkInv);
//check the resulting filters frequency response
efr(x)=firResp(wfk); //estimated freq. response
efrInv(x)=firResp(fkInv); //estimated freq. response length 2
plot(efr,efrInv,-.1,1.1,Pi+.1,-.5);
Band-pass/Band-reject out of a High-Pass and a Low-Pass
samplingF = 2048; //Hz
samplingTime = 1; //seconds
filterLength = 81; //points
shape(x) = if x < ( samplingF / 2 ) / 3 then 1 else 0; //say limit to a third of the interval
fr(x)=sampleN( shape, 0, 1, samplingF/2+1 ); //freq. response
//plot(fr);
fk(x)=ifft1(fr); //filter kernel 1
//plot(fk);
fkshr(x)=shrotS(fk,filterLength/2);
fktr(x)=truncateS(fkshr,0,filterLength-1);
fkn(x)=fktr(x)/sumS(fktr); //normalized kernel
blackmanW(x,filterLength) = 0.42 - 0.5 * cos( 2 * Pi * x / filterLength ) + 0.08 * cos( 4 * Pi * x / filterLength ); //Blackman
db(x)=8.685890*log(x);
wfk(x)=blackmanW(x/stepS(fkn),filterLength)*fkn(x); //the Blackman windowed filter
fkInv(x)=reverseFilter(wfk); //make high pass out of the low pass
plot(wfk,fkInv);
//check the resulting filters frequency response
efr(x)=firResp(wfk); //estimated freq. response
efrInv(x)=firResp(fkInv);
plot(efr,efrInv,-.1,1.1,Pi+.1,-.5);
//rf(x)=conv(wfk,fkInv); //band-pass -- THIS DOESN'T WORK, LOTS OF NOISE !!!
rf(x)=wfk(x)+fkInv(x); //result filter band-reject
rfn(x)=rf(x)/sumS(rf);
plot(rfn);
//rffr(x)=firResp(rfn);
//plot(rffr,-.1,1.1,Pi+.1,-.5);
//check the response but doing fft on the kernel
fknp(x)=padS(rf,2^nextpow2(samplingF*samplingTime)); //pad to the proper length
fkn1f(x)=fft1(fknp); //take fft
fkn1a(x)=abs(fkn1f(x))*samplingF/2; //take abs and normalize
fkn1db(x)=db(fkn1a(x)); //calculate in db
plot(fkn1a);
plot(fkn1db);
//Since the conv line above doesn't work we can get bandpass by inversion of bandreject
rfInv(x)=inverseFilter(rfn);
plot(rfInv);
rfInvFr(x)=firResp(rfInv);
plot(rfInvFr,-.1,1.1,Pi+.1,-.5);
Interactive test
samplingF = 2048; //Hz
samplingTime = 1; //seconds
controlLength = 81; //points
controlVar( controlLength, 50 );
nextOdd(x) = if isNextIntEven(x) then nextInt(x) + 1 else nextInt(x);
filterLength=nextOdd( controlLength );
rotLength=filterLength/2;
truncLength=filterLength-1;
cutoffF = .5; //fraction of the Niquist interval
controlVar(cutoffF,.5);
shape(x) = if x < ( samplingF / 2 ) * cutoffF then 1 else 0;
fr(x)=sampleN( shape, 0, 1, samplingF/2+1 ); //freq. response
//plot(fr);
fk(x)=ifft1(fr); //filter kernel 1
//plot(fk);
fkshr(x)=shrotS(fk,rotLength);
fktr(x)=truncateS(fkshr,0,truncLength);
fkn(x)=fktr(x)/sumS(fktr); //normalized kernel
blackmanW(x,filterLength) = 0.42 - 0.5 * cos( 2 * Pi * x / filterLength ) + 0.08 * cos( 4 * Pi * x / filterLength ); //Blackman
db(x)=8.685890*log(x);
wfk(x)=blackmanW(x/stepS(fkn),filterLength)*fkn(x); //the Blackman windowed filter
plot(wfk);
//check the resulting filters frequency response
efr(x)=firResp(wfk); //estimated freq. response
plot(efr,-.1,1.1,Pi+.1,-.5);
