Sign in
Search Online Books
Free Online Books
Chapters
Search Spectral Audio Signal Processing
Would you like to be notified by email when Julius Orion Smith III publishes a new entry into his blog?
Software Listing
The software listing below is in the SAIL programming language (generally a superset of ALGOL). Since procedures are defined before they are used, the listing starts out defining basic utilities, progressing to higher levels until finally the top level is reached. As a result, it is usually easier start read the top level procedure at the end first, and then work backwards to fill in the details.
The SAIL keyboard had various non-standard characters which are spelled out here using TeX equivalents such as &alpha#alpha; and &beta#beta;.
COMMENT Track instantaneous amplitudes and frequencies of multiple sinusoids; COMMENT Side-result: FFT-based filtering and/or inharmonic additive synthesis; COMMENT Completed in the summer of 1985 by Julius O. Smith III BEGIN "PARSHL" REQUIRE "{}<>" DELIMITERS; DEFINE #="comment",CrLf="'15&'12",tab={""&'11},Sam16="4",Cr="'15", ALT={'175&""},Thru={step 1 until}; REQUIRE "JOSLIB.REQ[LIB,JOS]" SOURCE!FILE; REQUIRE "REQUIRING UDP2:SIGLIB.REL[SIG,MUS] LIBRARY" MESSAGE; REQUIRE "UDP2:SIGLIB.REL[SIG,MUS]" LIBRARY; REQUIRE "JAMLIB.REL[SUB,SYS]" LIBRARY; REQUIRE "TRMTYP.SAI[LIB,JOS]" SOURCE!FILE; REQUIRE "DISPLA.REQ[LIB,JOS]" SOURCE!FILE; EXTERNAL PROCEDURE TRPINI(INTEGER CODE); # Enable floating point exception code; EXTERNAL PROCEDURE SUPCT; # Disable \alpha-I,\alpha-R,\alpha-T,\alpha-L etc; REQUIRE "RECORD.REQ[LIB,JOS]" SOURCE!FILE; # REQUIRE "MYIO.REQ[LIB,JOS]" SOURCE!FILE; COMMENT Declare File IO stuff; REQUIRE "FLTIO.REQ[LIB,JOS]" SOURCE!FILE; COMMENT Declare GetFlt; EXTERNAL BOOLEAN PROCEDURE AinInt(REFERENCE INTEGER Val; STRING Name); EXTERNAL BOOLEAN PROCEDURE AinReal(REFERENCE REAL Val; STRING Name); EXTERNAL STRING PROCEDURE Cvfs(REAL Val); EXTERNAL SIMPLE BOOLEAN PROCEDURE FNparse( STRING Arg; REFERENCE STRING Device; REFERENCE STRING Filename); EXTERNAL STRING PROCEDURE DEV(STRING Name); COMMENT Return DEVICE part of filename; EXTERNAL STRING PROCEDURE NAME(STRING Name); COMMENT Return NAME part of filename; EXTERNAL STRING PROCEDURE EXT(STRING Name); COMMENT Return extension part of filename; EXTERNAL STRING PROCEDURE PPN(STRING Name); COMMENT Return PPN part of filename; COMMENT SOUND IO DECLARATIONS; REQUIRE "SYS:DEFAUL.HDR" SOURCE!FILE; # REQUIRE "RHEAD.REL[SYS,MUS]" LOAD!MODULE; # REQUIRE "WHEAD.REL[SYS,MUS]" LOAD!MODULE; EXTERNAL BOOLEAN PROCEDURE SANDI( INTEGER Chan, StSamp, Nsamps; REFERENCE REAL ARRAY X; INTEGER FilSamps, Headed, FilPak, Xpack(-1)); EXTERNAL BOOLEAN PROCEDURE SANDO( INTEGER Chan, StSamp, Nsamps; REFERENCE REAL ARRAY X; REFERENCE INTEGER ARRAY Hist; REFERENCE INTEGER FilSamps; INTEGER Headed, FilPak, Xpack(-1)); # ReadH (rp, hdbuf, Fail, true if don't want printout); EXTERNAL PROCEDURE ReadH ( RECORD!POINTER (File) Rp; REFERENCE REAL ARRAY Hdbuf; REFERENCE INTEGER FAIL; BOOLEAN Silence (FALSE)); # WriteH (head, clock, pack, #chans, #ticks, maxamp, cstr); EXTERNAL PROCEDURE WriteH ( REFERENCE REAL ARRAY Head; REAL Clock; INTEGER Pack, #chans, #ticks; REAL Maxamp; STRING Cstr ); INTEGER ARRAY OutH,AmpH,FrqH[0:128]; REAL ARRAY Head[1:128]; INTEGER Nread,Brk; REAL Maxamp; BOOLEAN FAIL,Eof; STRING Ifile,Idev,Ofile,Odev; STRING AmpFile,AmpDev,FrqFile,FrqDev,Tstr; RECORD!POINTER (File) InFptr, OutFptr; RECORD!POINTER (File) AmpFptr, FrqFptr; DEFINE InF(x) = {File:x[InFptr]}, OutF(x) = {File:x[OutFptr]}; DEFINE AmpF(x) = {File:x[AmpFptr]}, FrqF(x) = {File:x[FrqFptr]}; COMMENT Filter input; # REQUIRE "JOSLIB.REQ[LIB,JOS]" SOURCE!FILE; # REQUIRE "RECORD.REQ[LIB,JOS]" SOURCE!FILE; # REQUIRE "FLTIO.REQ[LIB,JOS]" SOURCE!FILE; COMMENT Declare GetFlt; DEFINE MaxCoeffs = "8192"; COMMENT Maximum number of filter coefficients; STRING Ffile; INTEGER Ni,No; REAL ARRAY Ic,Oc[1:MaxCoeffs]; INTEGER Nfft,Nspec,Nh,Nx,Nhop,Ndec,i,Nhops; INTERNAL INTEGER Trace; DEFINE ShowSpectra="(Trace LAND 1)"; DEFINE TDpyFFA = " IF ShowSpectra THEN DpyFFA "; DEFINE ShowTime="(Trace LAND 2)"; DEFINE TDpyEd = " IF ShowTime THEN DpyEd "; DEFINE ShowPeakFinder="(Trace LAND 4)"; DEFINE Debug1="(Trace LAND 8)", Debug2="(Trace LAND 16)", Debug3="(Trace LAND 32)"; BOOLEAN HaveOfile,HaveIfile,HaveFfile,HavePack; STRING TmpStr; EXTERNAL INTEGER !SKIP!; IFC FALSE THENC # SAIL Esc-I interrupt facility; SIMPLE PROCEDURE EscI; # This is called upon <esc>I interrupt; START!CODE "EscI" TLZE '10,'400000; # If sign bit is on; MOVN '10,'10; # Convert sign-magnitude to 2' comp; MOVEM '10, Trace; # Save Ac 10; CALLI 0, '400024; # DISMIS (return to user level); END "EscI"; PROCEDURE Int!init; BEGIN EXTERNAL INTEGER JOBAPR; # JOB DATA AREA user interrupt routine pointer; JOBAPR <- LOCATION(EscI); ENABLE(15); # Enable user interrupt handling; Trace <- 0; END; COMMENT REQUIRE Int!Init INITIALIZATION; COMMENT rmoved by BIL because he thinks this is extremely dangerous -- you are depending on getting only esc-i interrupts and nothing else. The SAIL runtime world that tries to handle arithmetic exceptions (for example) depends on JOBAPR pointing to some procedure that can handle such interrupts (UDP2:TRIGS[NEW,SAI]/4P and UDP2:NWORLD[NEW,SAI]/31,32P); ELSEC IFC NOT DECLARATION(GOGTAB) THENC EXTERNAL INTEGER ARRAY GOGTAB[0:'300]; ENDC SIMPLE PROCEDURE Esci; Trace <- GOGTAB['200]; PROCEDURE Int!Init; BEGIN ENABLE(15); INTMAP(15,Esci,0); Trace <- 0; END; REQUIRE Int!Init INITIALIZATION; ENDC # Global declarations; BOOLEAN DoFlt,DoSynth,DoOut,SwapOut; # Flags for type of function wanted; REAL MinSep,Thresh,Hyst,DFmax1,DFmax2,DAmax,SigScl,Fc1,Fc2,Fs,DBscl; INTEGER MaxLins,MaxOscs,Nframes,Npartials,Frame,MinWid; BOOLEAN Quiet,InstantRise; COMMENT DerivedName - Generate name from root name and code string; STRING PROCEDURE DerivedName(STRING GivenName, CodeString); BEGIN "DerivedName" INTEGER i,j,ln,lc; STRING Gname,DevStr; STRING RootName; # REQUIRE "JOSLIB.REQ[LIB,JOS]" SOURCE!FILE; EXTERNAL STRING PROCEDURE DEV(STRING Name); COMMENT Return DEVICE part of filename; EXTERNAL STRING PROCEDURE NAME(STRING Name); COMMENT Return NAME part of filename; EXTERNAL STRING PROCEDURE EXT(STRING Name); COMMENT Return extension part of filename; EXTERNAL STRING PROCEDURE PPN(STRING Name); COMMENT Return PPN part of filename; Gname <- NAME(GivenName); ln <- LENGTH(Gname); lc <- LENGTH(CodeString); IF lc+ln LEQ 6 THEN RootName <- Gname&CodeString ELSE RootName <- Gname[1 to 6-lc]&CodeString; IF EQU(RootName,Gname) THEN BEGIN # Gak! Must generate a new name!; PRINT(" DerivedName got a name collision",CrLf); IF ln>1 THEN RootName <- Gname[1 to 6-lc-1]&CodeString&"2" ELSE RootName <- Gname[1 to 6-lc]&CodeString[1 to lc-1]&"2"; IF EQU(RootName,Gname) THEN RootName <- RootName[1 to 5]&"3"; END; DevStr <- Dev(GivenName); IF DevStr NEQ NULL THEN DevStr <- DevStr&":"; RETURN(DevStr&RootName&Ext(GivenName)&Ppn(GivenName)); END "DerivedName"; COMMENT Fprint - Print documentation from DSK; # Put this in JOSLIB; PROCEDURE Fprint(STRING Fname); COMMENT Print contents of file in fname; BEGIN STRING Help,Ttystr; INTEGER Fchan,Brk,Eof,Bt; BOOLEAN FAIL; DEFINE FF={('14&"")}; OPEN(Fchan <- GETCHAN,"DSK",0,2,0,2000,Brk,Eof); COMMENT Ascii input; LOOKUP(Fchan,Fname,FAIL); IF FAIL THEN USERERR(0,1," Can't find HELP file "&Fname); IF FAIL THEN RETURN; SETBREAK(Bt <- GETBREAK,FF,NULL,"sin"); DO BEGIN "pages" Help <- INPUT(Fchan,Bt); IF NOT EQU(Help[1 FOR 7],"COMMENT") THEN BEGIN MyOutDpy(Help,-3,-1); MyOutDpy ("Type RETURN to continue - ALT to return to command level",31,-1); MyOutDpy (Fname,33,-1); Ttystr <- INCHWL; Break!N; IF Ttystr="H" OR Ttystr = "?" THEN MyOutDpy(" This IS help!!",32,-1); IF !skip!="H" OR !skip! = "?" THEN MyOutDpy(" This IS help!!",32,-1); IF !skip! = ALT THEN BEGIN PRINT(CrLf); DONE "pages"; END; END; END "pages" UNTIL Eof; Relbreak(Bt); RELEASE(Fchan); END; COMMENT Help - Command Summaries; RECURSIVE PROCEDURE HELP(STRING Topic(NULL)); BEGIN "help" Break!N; IF EQU(Topic,"?") THEN BEGIN Fprint("PARSHL.JOS[MUS,DOC]"); PRINT(CrLf,CrLf," Type \alpha-? for command summary",CrLf); # PRINT(" Type <command!letter>\alpha-? for individual command description",CrLf); PRINT(" Type ?\alpha-? or ?,?<RETURN> for full documentation",CrLf,CrLf); RETURN; END; IF EQU(Topic,"TopLevel") THEN MyOutDpy(" Outer Level Command Summary [s> => string, [n> => integer, [d> => real number, \alpha- => CONTROL, and \beta- => META ALT - Exit PARSHL. ?\alpha-? - Print complete documentation. (Also ?,?[CR> does this.) [n>\alpha-W - Analysis Window (1=rect,2=triang,3:5=Hamming,6=Kaiser,7=Cheb) [n>\alpha-T - FFT length (must be a power of 2). [n>\alpha-D - Number of samples of input signal per FFT frame. [n>\alpha-H - Number of samples of input signal to advance between FFT's. [n>\alpha-L - Length of filter impulse response (if known). [n>\alpha-C - Decimation factor (1 => no decimation, 2 => every other sample, etc.) [s>\alpha-I - Input sound file. [s>\alpha-O - Output sound file. [n>\alpha-P - Output sound file packing mode. (\alpha-P with no arg explains codes.) [s>\alpha-F - Input filter file. ; - Comment. [ESC>[n>I codes: (n can be any sum of the following): [n> display(s) enabled --- ------------------ 0 none (type [ESC>I to turn off all running displays) 1 time buffers ([ESC>1I) 2 spectral buffers ([ESC>2I) 4 peak-finder 8 debug level 1 16 debug level 2 32 debug level 3 ------------------------------------------------------------------------ ",-3,-1); IF EQU(Topic,"SubLevel") THEN MyOutDpy(" Inner Level Command Summary [s] => string, [n] => integer, [d] => real number, \alpha- => CONTROL, and \beta- => META ALT - Exit to Outer level. ?\alpha-? - Print complete documentation. (Also ?,?[CR] does this.) [d]\alpha-M - Minimum spacing between partials in Hz [n]\alpha-N - Maximum number of partials to track [n]\alpha-O - Maximum number of oscillators to use [d]\alpha-F - First (lowest) partial frequency allowed (Hz) [d]\alpha-L - Last (highest) partial frequency allowed (Hz) [d]\alpha-T - Peak detection threshold in dB [d]\alpha-H - Spectral wiggle tolerance in dB [d]\alpha-D - Maximum allowed frequency deviation per frame in Hz at LOW frequency [d]\alpha-U - Maximum allowed frequency deviation per frame in Hz at HIGH frequency [n]\alpha-B - Buffer display code. If negative, all time buffers are displayed. \alpha-S - Show spectra of filter input and output on each FFT hop. (\beta-S disables.) ; - Comment. ------------------------------------------------------------------------ ",-3,-1); END "help"; COMMENT DpyFFA - Display interleaved dB spectra as computed by FFA in SIGLIB; PROCEDURE DpyFFA(REAL ARRAY S; INTEGER Nfft; STRING T; REAL Fs(1)); COMMENT For displaying interleaved dB spectra as computed by FFA in SIGLIB; BEGIN "DpyFFA" INTEGER i,Nspec; REAL ARRAY Spec[1:(Nfft DIV 2)+1]; Nspec <- (Nfft DIV 2) + 1; FOR i <- 1 STEP 1 UNTIL Nspec DO Spec[i] <- (10.0/LOG(10.0))*LOG( (S[2*i-1]^2+S[2*i]^2) MAX 1.0@-35); DpyEd(Spec,Nspec,T,"MAGNITUDE (DB)",0,Fs/2); END "DpyFFA"; COMMENT Qinterp - Quadratic interpolation; INTERNAL SIMPLE REAL PROCEDURE Qinterp(REAL Ym1,Y0,Yp1; BOOLEAN InteriorX(TRUE)) ; COMMENT Fit a parabola Y[X] = A*X^2+B, through the three equally spaced Y-values Ym1 = Y[-1], Y0 = Y[0], and Yp1=Y[1], and return the X-value where the extremum is attained. For example, if -1 is returned, then Y[-1] is the point of zero slope in the parabola through the three points. If InteriorX is TRUE then if the extremum lies outside the interval [-1,1], a warning is printed and the returned value is clipped to lie at either 1 or -1 as appropriate. ; BEGIN "Qinterp" REAL X; X <- (Yp1 - Ym1)/(2*(2*Y0 - Yp1 - Ym1)); IF InteriorX AND (ABS(X)>1) THEN PRINT(" Qinterp: Clipping analytic extremum to ", X <- (IF X>0 THEN 1 ELSE -1),CrLf); IF Debug3 THEN PRINT("Qinterp: Given Y's ",Ym1," ",Y0," ",Yp1, ", extremum is at X = ",X,CRLF); RETURN(X); END "Qinterp"; COMMENT FindPeaks - Generic peak finder; IFC NOT DECLARATION(IntNoSet) THENC DEFINE IntNoSet={(1 LSH 35)}; IFC NOT DECLARATION(RealNoSet) THENC DEFINE RealNoSet={(-1.0@35)}; INTEGER PROCEDURE FindPeaks(REAL ARRAY X,Peaks,PeakLocs; REAL Thresh(RealNoSet),Hyst(RealNoSet); INTEGER MaxPeaks(IntNoSet),MinWidth(IntNoSet),BlankWidth(IntNoSet), I1(IntNoSet),I2(IntNoSet)); COMMENT Find amplitudes and indices of all peaks in X = array of values to be searched for peaks. (It is overwritten.) Peaks[1:MaxPeaks] = peak amplitudes. PeakLocs[1:MaxPeaks] = peak indices. Thresh = Threshold below which no peaks are to be sought Hyst = Wiggles less than Hyst in magnitude are ignored. MaxPeaks = Maximum number of peaks to be found (starting w largest) MinWidth = Minimum peak width in indices, after blanking. Narrower peaks are removed and ignored. BlankWidth = Minimum blanking interval in indeces. If = IntNoSet, then blank to nearest local minimum to within Hyst interval. Blanking is done before width is measured for efficiency reasons (would rather measure width first). i1,i2 = index boundaries for the search within X. The return value is the number of peaks actually found. ; BEGIN "FindPeaks" INTEGER i,j,Npeaks,NdxReach,Owid,Odig,Peak,Poff,PLoff,M1,M2,Nfound,Lb,Ub; REAL Xmax,Xmin; GETFORMAT(Owid,Odig); SETFORMAT(0,5); NdxReach <- (IF (BlankWidth NEQ IntNoSet) THEN (BlankWidth-1)/2 MAX 1 ELSE 0); Poff <- ARRINFO(Peaks,1)-1; PLoff <- ARRINFO(PeakLocs,1)-1; Npeaks <- (IF MaxPeaks NEQ IntNoSet THEN MaxPeaks ELSE ARRINFO(Peaks,2)-Poff); Npeaks <- Npeaks MIN (ARRINFO(Peaks,2)-Poff) MIN (ARRINFO(PeakLocs,2)-PLoff); Lb <- ARRINFO(X,1); Ub <- ARRINFO(X,2); IF I2=IntNoSet THEN I2 <- ARRINFO(X,2); IF I1=IntNoSet THEN I1 <- ARRINFO(X,1); I1 <- (Lb MAX I1 MIN Ub); I2 <- (I1 MAX I2 MIN Ub); ArrMin(Xmin,X,I1,I2); ArrMax(Xmax,X,I1,I2); IF Thresh=RealNoSet THEN Thresh <- Xmin; IF Hyst=RealNoSet THEN Hyst <- (Xmax-Xmin)/100; M1 <- M2 <- 0; COMMENT [m1,m2] = index interval of last peak; Nfound <- 0; FOR Peak <- 1 STEP 1 UNTIL Npeaks DO BEGIN "fp" INTEGER MaxLoc,TmpI; REAL MaxVal,ClobVal,TmpR; IF M1=I1 AND M2=I2 THEN DONE "fp"; MaxLoc <- ArrMax(MaxVal,X,I1,I2); IF MaxVal<Thresh THEN DONE "fp"; Nfound <- Nfound + 1; PeakLocs[Nfound+PLoff] <- MaxLoc + Qinterp(X[(MaxLoc-1) MAX I1],MaxVal,X[(MaxLoc+1) MIN I2]); IF MaxLoc=I1 AND Debug1 THEN PRINT("*** Peak is at right of search interval for peak ",Nfound," ***",CrLf); IF MaxLoc=I2 AND Debug1 THEN PRINT("*** Peak is at left of search interval for peak ",Nfound," ***",CrLf); Peaks[Nfound+Poff] <- MaxVal; COMMENT Now slice off peak so we don't find it again; M1 <- (MaxLoc-NdxReach) MAX I1; M2 <- (MaxLoc+NdxReach) MIN I2; ArrMin(ClobVal,X,M1,M2); TmpR <- X[M1]; WHILE M1>I1 AND TmpR+Hyst GEQ X[M1-1] DO BEGIN TmpR <- TmpR MIN X[M1-1]; M1 <- M1-1 END; ClobVal <- ClobVal MIN TmpR; TmpR <- X[M2]; WHILE M2<I2 AND TmpR+Hyst GEQ X[M2+1] DO BEGIN TmpR <- TmpR MIN X[M2+1]; M2 <- M2+1 END; ClobVal <- ClobVal MIN TmpR; FOR i <- M1 STEP 1 UNTIL M2 DO X[i] <- ClobVal; IF (M2-M1+1 < MinWidth) OR MaxLoc=I1 OR MaxLoc=I2 # JOS/10/25/85; THEN BEGIN "Throw it back" Nfound <- Nfound-1; CONTINUE "fp"; END "Throw it back"; IF ShowPeakFinder THEN DpyEd(X,I2-ARRINFO(X,1)+1,"X remains after peak "& CVS(Nfound)&" at i="&CVS(MaxLoc)); END "fp"; SETFORMAT(Owid,Odig); RETURN(Nfound); END "FindPeaks"; COMMENT FindPartials - Find all significant spectral lines for current frame; PROCEDURE FindPartials(REFERENCE INTEGER Npartials; REAL ARRAY XmagDB,LinAmp,LinFrq; REAL Fs,MinSep,Thresh,Hyst; INTEGER Nfft,Frame,MinWid; REAL Fc1(RealNoSet),Fc2(RealNoSet)); COMMENT Find amplitudes (dB) and frequencies (Hz) of all partials in XmagDB. Npartials = the maximum number of partials to find. It is set to the actual number of partials found on return. XmagDB[1:Nfft/2+1] = DB magnitude spectrum LinAmp[1:Nharms,1:Nframes] = Partial amplitudes LinFrq[1:Nharms,1:Nframes] = Partial frequencies Fs = sampling rate (Hz) MinSep = minimum line spacing (Hz) Thresh,Hyst are args to FindPeaks Nfft = Size of FFT used in computing XmagDB Frame = current time frame number MinWid = minimum acceptable peak width in FFT bins. Fc1,Fc2 = lower and upper frequency bounds in Hz. No partials will be allowed outside this interval. ; BEGIN "FindPartials" REQUIRE "HACKS.REL[SUB,SYS]" LIBRARY; EXTERNAL SIMPLE PROCEDURE QISort(REAL ARRAY A,Data; INTEGER LB,UB); INTEGER i,j,Nspec,BinInt,BinReach,Owid,Odig,Partial,MinWid; REAL ARRAY Xsave[1:Nspec <- Nfft/2+1]; GETFORMAT(Owid,Odig); SETFORMAT(0,7); BinInt <- (MinSep/Fs)*Nfft+.5; # Blanking interval for spectral peak (in bins); IF Frame=1 THEN PRINT(CrLf," Blanking interval is ",BinInt," bins",CrLf); ARRTRAN(Xsave,XmagDB); # Save input array; Npartials <- FindPeaks(XmagDB,LinAmp,LinFrq,Thresh,Hyst, Npartials,MinWid,BinInt,((1+Nfft*Fc1/Fs) MAX 1),((1+Nfft*Fc2/Fs) MIN Nspec)); ARRTRAN(XmagDB,Xsave); # Sort for ascending frequency (for convenience only); IF Npartials>0 THEN QIsort(LinFrq,LinAmp,1,Npartials); IF Debug2 THEN PRINT(CrLf); FOR Partial <- 1 STEP 1 UNTIL Npartials DO BEGIN "fp" LinFrq[Partial] <- Fs*(LinFrq[Partial]-1)/Nfft; IF Debug2 GEQ Partial THEN PRINT("Frame=",Frame," Amp(dB)=",LinAmp[Partial], " Freq(Hz)=",LinFrq[Partial],"=",Partial,"*",MinSep,"*", LinFrq[Partial]/(Partial*MinSep),CrLf); END "fp"; IF ShowPeakFinder THEN BEGIN "ShowHarms" REAL ARRAY CFrqs[1:Nspec]; REAL Val; Val <- LinAmp[1]; j <- 1; FOR Partial <- 1 STEP 1 UNTIL Npartials DO BEGIN FOR i <- j STEP 1 UNTIL Nfft*(LinFrq[Partial]/Fs) DO CFrqs[i] <- Val; j <- i; Val <- LinAmp[Partial]; END; Dpy2(XmagDB,CFrqs,i, "Centers of pitch-wide search regions for frame "&CVS(Frame), "MAGNITUDE (DB)",0,Fs/2,-30,30); END "ShowHarms"; SETFORMAT(Owid,Odig); END "FindPartials"; COMMENT OutPartials - Write partial amplitudes and frequencies to disk; PROCEDURE OutPartials(INTEGER Npartials,Nframes; REAL ARRAY Amps,Frqs; REAL Fs,Thresh); BEGIN "OutPartials" INTEGER N,AmpP,FrqP; N <- Npartials*Nframes; Sando(AmpF(Channel),AmpP,N,Amps,AmpH,AmpP,TRUE,AmpF(Pack)); Sando(FrqF(Channel),FrqP,N,Frqs,FrqH,FrqP,TRUE,FrqF(Pack)); Sando(AmpF(Channel),0,0,Amps,AmpH,AmpP,TRUE,AmpF(Pack)); Sando(FrqF(Channel),0,0,Frqs,FrqH,FrqP,TRUE,FrqF(Pack)); END "OutPartials"; COMMENT Decimate - Downsample by integer factors; INTEGER PROCEDURE Decimate(REAL ARRAY A; INTEGER N,M,I(0)); COMMENT Downsample array A[1:N] by M. I nonzero means initialize; COMMENT Return the number of samples produced; BEGIN "Decimate" OWN INTEGER P; INTEGER j; IF M LEQ 1 THEN RETURN(N); IF I NEQ 0 OR P LEQ 0 THEN P <- 1; j <- 0; FOR i <- P STEP M UNTIL N DO A[j <- j+1] <- A[i]; P <- i-N; RETURN(j); END "Decimate"; COMMENT GetWin - Compute overlap-add analysis window; INTEGER WinType; DEFINE Rectangular="1", Triangular="2", Hamming="3", GenHamming="4", Hanning="5", Kaiser="6", Chebyshev = "7", Nwins = "7"; PRELOAD!WITH "Rectangular", "Triangular", "Hamming", "GenHamming", "Hanning", "Kaiser", "Chebyshev"; STRING ARRAY WinNames[1:Nwins]; DEFINE WinStr = "WinNames[WinType]"; PROCEDURE GetWin(REAL ARRAY W; INTEGER Wtype,Nw; REAL P3(-1.),P4(-1.)); COMMENT Generate analysis window taking special care to ensure that it will overlap-add to unity in the case of a Hamming window with hopsize = Nw/(2K) Wtype Window ----- ------ 1 Rectangular 2 Triangular 3 Hamming 4 Generalized Hamming 5 Hanning 6 Kaiser 7 Chebyshev ; BEGIN "GetWin" # REQUIRE "SIGLIB.REL[SUB,SYS]" LIBRARY; EXTERNAL PROCEDURE !WINFLT(REAL ARRAY H; REFERENCE INTEGER NH; REFERENCE INTEGER WINTYP; REFERENCE INTEGER BNDTYP; REAL ARRAY P); REAL ARRAY WinPars[1:4]; WinPars[1] <- 0; # Lower cutoff frequency; WinPars[2] <- 0; # Upper cutoff frequency; IF WinType=Kaiser AND P3<0 THEN DO BEGIN IF NOT AinReal(P3 <- 60,"Kaiser stopband rejection in DB") THEN RETURN; IF P3 < 0 THEN PRINT("Stopband rejection must be POSITIVE dB...like 60",CRLF); END UNTIL P3 GEQ 0; IF WinType=GenHamming AND P3<0 THEN DO IF NOT AinReal(P3 <- 0.54,"Alpha (0:4)") THEN RETURN UNTIL (0 LEQ P3 LEQ 4); IF WinType=Chebyshev AND P4 LEQ 0 THEN DO BEGIN "GetCheb" IF NOT AinReal(P3, "Chebyshev stopband rejection in DB (or minus transition width in Hz/Srate)") THEN RETURN; IF P3<0 THEN BEGIN P4 <- P3; P3 <- 0; END ELSE P4 <- 0; END "GetCheb" UNTIL (0 LEQ P4<1/2 AND 0 LEQ P3); WinPars[3] <- P3; WinPars[4] <- P4; IF (Nw MOD 2) = 1 THEN BEGIN REAL ARRAY TmpBuf[1:2*Nw+1]; !WinFlt(TmpBuf,Nw,Wtype,1,WinPars); # Analysis window = lowpass, Fc=0; ARRTRAN(W,TmpBuf); IF Wtype>Rectangular THEN W[Nw] <- 0; # For odd lengths, last sample zeroed; END ELSE BEGIN REAL ARRAY TmpBuf[1:2*Nw+1]; INTEGER i; !WinFlt(TmpBuf,i <- 2*Nw+1,Wtype,1,WinPars); FOR i <- 1 STEP 1 UNTIL Nw DO W[i] <- TmpBuf[2*i]; END; BEGIN "nrmlz" # REQUIRE "JOSLIB.REQ[LIB,JOS]" SOURCE!FILE; EXTERNAL REAL PROCEDURE MaxArr(INTEGER n; REAL ARRAY y); INTEGER i; REAL Wmax,Wscl;