Digital Filter Sound Effects

Donald Daniel

July 2017 revised mar 2020

up one level

If the lines of text are too long you can fix the problem with these instructions.

Below is a program that is an example of creating sound effects with digital filters.

The last time I heard a B36 airplane was in about 1957. I am writing this in 2017, so that was 60 years ago. The recordings I have heard do not have as much low frequency as I remember. So I wrote a program to simulate the sound I remember. I did not succeed as well as I had hoped. Perhaps the technique will be of interest to someone else.

Suppose the sound you want consists of a thump, a chirp and a rattle all at the same time. I think this technique will require a different filter for each component of a complicated sound. Then the sounds generated by each filter should be added in the proper proportions to achieve the desired sound. But that is not what I did here. I only used one filter.

The program sends a stream of impulses through a digital filter. The repetition rate of the stream of impulses is 36Hz, about the repetition rate of the propeller blades at cruise rpm. The program allows the user to enter poles and zeroes. First an impulse is sent through the filter once to find the amplitude of the response. This is used to normalize the gain. Then a stream of impulses is sent through the filter to get the desired sound waveform. Then it is sent through a highpass filter to eliminate the DC component which is not present in any real sound wave. It is computed in floating point, then converted to 16 bit integer for conversion to wav format.

When the program is run, enter "cp" for a complex pole pair. Then -140 1400. Then "q3" to start the calculation.

The program uses separately compiled modules listed in the article "example digital filters" at this website.

To download the file temp.wav click here.

up one level

MODULE b36a;
IMPORT Msg,In,Out,Files,BinaryRider,rm:=RealMath,
cx:=cxarith,cmpt, OS:ProcessManagement;

CONST pi=3.14159;rate=48.0E3;strlngth=7;
q=1;er=2;rp=3;rz=4;zo=5;cp=6;cz=7;q3=8;
last=9;freq=45;fhp=10.0;

TYPE
str=ARRAY strlngth OF CHAR;
cma=ARRAY last+1 OF str;
 
VAR
t,max,est,xkm1,ykm1,gf:REAL;first:BOOLEAN;
cmdara:cma;command,samples,g:LONGINT;
outfile:BinaryRider.Writer;
resv:Msg.Msg; outvar:Files.File;fr:cmpt.fltrrec;

PROCEDURE initvar;
BEGIN
cmdara[q      ]:='q';
cmdara[er      ]:='er';
cmdara[rp      ]:='rp';
cmdara[rz      ]:='rz';
cmdara[cp      ]:='cp';
cmdara[cz      ]:='cz';
cmdara[zo      ]:='zo';
cmdara[cz      ]:='cz';
cmdara[q3      ]:='q3';
cmdara[last    ]:='last';
first:=TRUE;t:=1.0/rate;
samples:=ENTIER(rate/freq);
outvar:=Files.New('temp.s16',{Files.write},resv);
outfile:=BinaryRider.ConnectWriter(outvar);
END initvar;


PROCEDURE readwd(VAR c:str);
BEGIN
In.Identifier(c);
IF In.Done() & (c # "") THEN
ELSE HALT(1) END;
END readwd;

PROCEDURE determine(VAR cmdvar:LONGINT);
CONST bell=7;
VAR cmi:LONGINT;strvar:str;
BEGIN
Out.String('enter command');Out.Ln;
cmdvar:=er;
readwd(strvar);
FOR cmi:=q TO last DO
IF (cmdara[cmi]=strvar)THEN cmdvar:=cmi;END;END;
IF cmdvar=er THEN
Out.Char(CHR(bell));Out.String('<--<< error');Out.Ln;END;
END determine;

PROCEDURE setf1proc;
BEGIN
Out.String('enter f1 normalization frequency');Out.Ln;
Out.String('at a freq with nonzero response');Out.Ln;
In.Real(fr.f1);
fr.f1:=2*pi*fr.f1;
END setf1proc;

PROCEDURE newpz;
VAR i,j:LONGINT;
BEGIN
i:=0;j:=0;fr.nzo:=0;
REPEAT
Out.String('enter q3 to quit this menu');Out.Ln;
Out.String('enter choice:');Out.Ln;
Out.String('rp rz zo cp cz');Out.Ln;
determine(command);
CASE command OF
er:;
|rp:Out.String('enter negative location');Out.Ln;
INC(i); In.Real(fr.p[i].sigma);fr.p[i].omega:=0.0;
|rz:Out.String('enter location');Out.Ln;
INC(j); In.Real(fr.z[j].sigma);fr.z[j].omega:=0.0;
|zo:setf1proc;Out.String('zero at origin created');Out.Ln;
INC(fr.nzo);
INC(j); fr.z[j].sigma:=0.0;fr.z[j].omega:=0.0;
|cp:Out.String('enter real then imaginary parts');Out.Ln;
INC(i); In.Real(fr.p[i].sigma);In.Real(fr.p[i].omega);
INC(i); fr.p[i].sigma:=fr.p[i-1].sigma;
fr.p[i].omega:=-fr.p[i-1].omega;
|cz:Out.String('enter real then imaginary parts');Out.Ln;
INC(j); In.Real(fr.z[j].sigma);In.Real(fr.z[j].omega);
INC(j); fr.z[j].sigma:=fr.z[j-1].sigma;
fr.z[j].omega:=-fr.z[j-1].omega;
|q3:;
END(*case*);
UNTIL command=q3;fr.n:=i;fr.nz:=j;
END newpz;

PROCEDURE normalize;
VAR i,k:LONGINT;xk,yk:cx.complex;
BEGIN
max:=0.0;
FOR k:=1 TO fr.n DO fr.p[k].first:=TRUE;
fr.p[k].ykm1:=cx.zero;END;
FOR k:=1 TO fr.nz DO fr.z[k].first:=TRUE;
fr.z[k].xkm1:=cx.zero;END;
FOR i:=1 TO samples DO
IF i=1 THEN xk:=cx.one ELSE xk:=cx.zero END;
cmpt.pzfltr(t,xk,yk,fr); 
IF cx.abs(yk)>max THEN max:=cx.abs(yk);END;END;
END normalize;

PROCEDURE highpass(t,fhp,xk:REAL;VAR yk:REAL);
(*Variables only used inside procedure must be declared
globally: first,est,gf,xkm1,ykm1. first must be initialized
globally. With filter 3dB point fhp at 10Hz, the result
of the filter is that the DC component reduced to a tenth,
20dB down, at 0.04 seconds*)
BEGIN
IF first THEN est:=rm.exp(-2*pi*fhp*t);
gf:=(1+est)/(1-est);first:=FALSE; 
xkm1:=0.0;ykm1:=0.0; END(*IF*);
yk:=(xk-xkm1)/2;xkm1:=xk; xk:=yk;
yk:=est*ykm1+(1-est)*xk;ykm1:=yk; yk:=gf*yk;
END highpass;

PROCEDURE gen;
VAR i,k:LONGINT;yki:INTEGER;
xkr,ykr:REAL;xk,yk,mplscx:cx.complex;
BEGIN
cx.rmult(1/max,cx.one,mplscx);
FOR k:=1 TO fr.n DO fr.p[k].first:=TRUE;
fr.p[k].ykm1:=cx.zero;END;
FOR k:=1 TO fr.nz DO fr.z[k].first:=TRUE;
fr.z[k].xkm1:=cx.zero;END;
FOR i:=0 TO ENTIER(3*rate) DO
IF (i MOD samples)=0 THEN xk:=mplscx ELSE xk:=cx.zero END;
cmpt.pzfltr(t,xk,yk,fr);xkr:=yk.r; 
highpass(t,fhp,xkr,ykr);
yki:=SHORT(ENTIER(16.0E3*ykr));
IF i>2000 THEN outfile.WriteInt(yki);END;
END(*i*);
(*outvar.Close;*)
END gen;

BEGIN
initvar;newpz;normalize; gen;
g:=ProcessManagement.system("sox -r 48k temp.s16 temp.wav");
g:=ProcessManagement.system("play temp.wav");
outvar.Close;
END b36a.


up one level

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.