Digital Filter Sound Effects

Donald Daniel

July 2017

up one level

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

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 -2000 2000. Then "cp" for another pole pair. Then -41 1111. 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 b36.wav click here.

up one level

MODULE b36a;
(*written 2017 by donald daniel.  
This program is free software:
you can redistribute it and/or modify it under the terms of the
GNU General Public License as published by the Free Software
Foundation, version 3 of the License.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General
Public License along with this program.  If not, see
http://www.gnu.org/licenses/.  *)

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=36;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('b36.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 b36.s16 b36.wav");
g:=ProcessManagement.system("play b36.wav");
g:=ProcessManagement.system("aac-enc b36.wav b36.aac");
END b36a.

up one level