Below is the original source listing for resample. It is written in the SAIL programming language, with time-critical inner loops written in FAIL (the assembly language used for the PDP KL-10 at the Stanford AI Lab at the time).
This source code is released under the BSD License
Note that SAIL's left-arrow has been replaced by underbar (_) below.
BEGIN "SRC" COMMENT Implement sampling rate conversions by (almost) arbitrary factors. Based on SRCONV.SAI[SYS,MUS] with new algorithm by JOS&PG. The program internally uses 16-bit data and 16-bit filter coefficients. Reference: "A Flexible Sampling-Rate Conversion Algorithm," J. O. Smith and P. Gossett, ICASSP, San Diego, 1984. To do: Fast 16-bit IO not yet implemented (buggy routines have been written but are not enabled). Also, need to add stereo case. Fast 16-bit stereo IO routines by MMM are here but not yet used. * Need overflow detection in fail section. Also want to saturate instead of wrap-around and report number of clipped samples. Modification history: 8/27/83 - JOS - Changed filter file format (added Nwing). 8/29/83 - JOS - Removed 0.95 guard factor in filter normalization. 8/30/83 - JOS - Changed width of scale factor (NlpScl) from 8 to 17. Moved conversion constants to SRC.PRM[SYS,MUS]. 9/ 3/83 - JOS - Changed ASH to LSH in filter IO to get 1 more MSB. 2/11/85 - JOS - Absorbed REQUIRE file "SRC.PRM" (page 10). Only JOS's XF software uses it besides SRC (as far as I know). 5/ 2/85 - JOS - Changed Xoff from ((Nmult-1)/2)*... to ((Nmult+1)/2)*.... ; REQUIRE "{}<>" DELIMITERS; DEFINE #="comment",CrLf="'15&'12",tab={""&'11},Sam16="4",Cr="'15", TwoTracks={(2*31*128)}, thru={step 1 until}; DEFINE Saf = {}; REQUIRE "SYS:DEFAUL.HDR" SOURCE!FILE; REQUIRE "JOSLIB.REQ[LIB,JOS]" SOURCE!FILE; REQUIRE "JAMLIB.REL[SUB,SYS]" LIBRARY; 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)); EXTERNAL PROCEDURE FIRLPF(REAL ARRAY Imp; INTEGER N; REAL Cutoff,Clock,Beta); EXTERNAL STRING PROCEDURE DEFAULT(STRING s; REFERENCE STRING Dev; STRING Defext,Defdev); EXTERNAL PROCEDURE TRPINI(INTEGER Flag); INTEGER PROCEDURE Pwr2(INTEGER N); BEGIN "Pwr2" INTEGER Nfft; Nfft _ LOG(N)/LOG(2) - 0.000001; Nfft _ 2^(Nfft+1); RETURN(Nfft); END "Pwr2"; # Display lowpass filter; INTERNAL INTEGER TRACE; PROCEDURE CkFIR(INTEGER ARRAY Imp; INTEGER Nwing; REAL Dh; BOOLEAN Interp); COMMENT Display frequency response of interpolation lowpass; BEGIN "CkFIR" # REQUIRE "JOSLIB.REQ[LIB,JOS]" SOURCE!FILE; EXTERNAL PROCEDURE FFT(INTEGER Nfft; REAL ARRAY X,A,B); EXTERNAL PROCEDURE ShowLogMag(INTEGER Nfft; REAL ARRAY A,B; STRING Title); REDEFINE NoSet="(1.7014118@38)"; EXTERNAL INTEGER PROCEDURE DPYED (REAL ARRAY Data; INTEGER Npoints(-1); STRING Xlabel(NULL),Ylabel(NULL); REAL Vxmin(NoSet),Vxmax(NoSet), Miny(NoSet),Maxy(NoSet)); PROCEDURE FFTwing(INTEGER Nwing,N; REAL Saf ARRAY X); COMMENT Given right wing of FIR filter in X[1:Nwing], create left wing at right of X[1:N] suitable for FFT (for frequency response); BEGIN "FFTwing" INTEGER i; REAL Saf ARRAY Temp[1:Nwing]; ARRTRAN(Temp,X); # Temp _ right wing; FOR i_1 Thru Nwing-1 DO X[N-i+1]_X[i+1]; # Far right _ left wing; FOR i_1 Thru Nwing DO X[i]_Temp[i]; # Far left _ right wing; FOR i_Nwing+1 Thru N-Nwing+1 DO X[i]_0; # Zero middle; END "FFTwing"; INTEGER Nfft,Nspec; Nfft _ Pwr2(2*Nwing/Dh)*4; Nspec _ (Nfft DIV 2) + 1; BEGIN INTEGER k,Iho; REAL Ho; REAL ARRAY X[1:Nfft],A,B[1:Nspec]; k_0; IF Interp THEN FOR Ho_0 STEP Dh UNTIL Nwing-1 DO BEGIN Iho_Ho; X[k_k+1]_Imp[Iho]+(Ho-Iho)*(Imp[Iho+1 MIN Nwing-1]-Imp[Iho]); END ELSE FOR Ho_0 STEP Dh UNTIL Nwing-1 DO X[k_k+1]_Imp[Ho]; DpyEd(X,k,"Right wing of effective impulse response subsampled by "&CVF(Dh)); FFTwing(k,Nfft,X); FFT(Nfft,X,A,B); ShowLogMag(Nfft,A,B,"Frequency response "&(IF Interp THEN "WITH" ELSE "WITHOUT")&" interpolation"); END; END "CkFIR"; # Two-channel 16-bit IO (by MMM) (not yet installed) These take the place of SANDI and SANDO. The use of buffered output without any USETOs allows writing to the tape drive, and the reading of samples and splitting into 2 arrays is more efficiently accomplished. ; PROCEDURE I16B2C (INTEGER Chan, StartSamp, Nsamps; Saf INTEGER ARRAY Chan1, Chan2); BEGIN "I16B2C" INTEGER I, J, V1, V2, Rec; Saf INTEGER ARRAY InArr [1:Nsamps]; StartSamp _ StartSamp + 128; Rec _ (StartSamp / 128); USETI (Chan, Rec+1); J _ StartSamp - (Rec*128); ARRYIN (Chan, InArr[1], J); ARRYIN (Chan, InArr[1], Nsamps); FOR I _ 1 STEP 1 UNTIL Nsamps DO BEGIN J _ InArr[I]; Chan1[I] _ IF (J LAND '20000000000) THEN (J LSH -16) LOR '777777600000 ELSE (J LSH -16) LAND '177777; Chan2[I] _ IF (J LAND '100000) THEN J LOR '777777600000 ELSE J LAND '177777; END; END "I16B2C"; PROCEDURE W16B2C (INTEGER Chan, Nsamps; Saf INTEGER ARRAY Chan1, Chan2); BEGIN "W16B2C" INTEGER I, A1, A2; Saf INTEGER ARRAY OutArr [1:Nsamps]; FOR I _ 1 STEP 1 UNTIL Nsamps DO OutArr[I] _ (Chan1[I] LSH 16) LOR (Chan2[I] LAND '177777); ARRYOUT(Chan, OutArr[1], Nsamps); END "W16B2C"; # Single channel 16-bit IO. I16B1C seems to work only on the first call to it. W16B1C has to be rewritten because it only writes multiples of 128 words. ; PROCEDURE I16B1C (INTEGER Chan, StartSamp, Nsamps; Saf INTEGER ARRAY Data); BEGIN "I16B1C" INTEGER I, J, K, Rec,Ns,Ss,Nz; Saf INTEGER ARRAY InArr [1:(Nsamps DIV 2) + 1 + 128]; Ns _ (Nsamps DIV 2) + 1 + 128; # Enough for partial record read; IF StartSamp<0 THEN BEGIN Nz_-StartSamp; StartSamp_0; END ELSE Nz _ 0; Ss _ StartSamp DIV 2; # True word number containing first sample; Ss _ Ss + 128; # Step over header; Rec _ (Ss DIV 128); # First Record number, numbering from 0; USETI (Chan, Rec+1); # System numbers from 1; J _ Ss - (Rec LSH 7); # Ss - Rec*128 = No. leading words to ignore; ARRYIN (Chan, InArr[1], Ns); PRINT(" May need to save and restore AC's here",CrLf); QUICK!CODE "In16c1" DEFINE A="1",D="2",De="3",T="4"; LABEL LH,RH,DONE; Move D,Data; # Pointer to current Data sample; Add D,Nz; # Number of initial zeros; Move De,Data; # Pointer to end of Data + 1; Add De,Nsamps; # An absolute address; Move A,InArr; # Pointer to current input array sample; Add A,J; # Step over extra samples; Move T,StartSamp; # 0 is first sample of file; Trne T,1; # Odd?; Jrst RH; # Yes, first word is in right half; LH: Caml D,De; # Data filled?; Jrst DONE; # Yes, return; Move T,@A; # Get current input word; LSH T,4; # Left justify it so sign bit is in bit 0; ASH T,-20; # Right justify with sign extension; Movem T,@D; # Deposit into Data array; Addi D,1; # and advance output pointer; RH: Caml D,De; # Data filled?; Jrst DONE; # Yes, return; Move T,@A; # Get current input word; LSH T,20; # Left justify it so sign bit is in bit 0; ASH T,-20; # Right justify with sign extension; Movem T,@D; # Deposit into Data array; Addi D,1; # Advance output pointer; Aoja A,LH; # Advance input pointer; DONE: Jfcl; # Nothing left to do I hope; END "In16c1"; END "I16B1C"; PROCEDURE W16B1C (INTEGER Chan, Nsamps; Saf INTEGER ARRAY Data); BEGIN "W16B1C" INTEGER I, J, Nout; OWN INTEGER Last; OWN BOOLEAN LastX; Saf INTEGER ARRAY OutArr [1:Nsamps]; QUICK!CODE "O16c1" DEFINE A="1",D="2",De="3",T="4"; LABEL LH,RH,LDone,Rdone,DONE; Move D,Data; # Pointer to current Data sample; Move De,D; # Pointer to end of Data + 1; Add De,Nsamps; # An absolute address; Move A,OutArr; # Pointer to current output array sample; Move T,Nsamps; # 0 is first sample of file; Skipn LastX; # Is there a leftover from last call?; Jrst LH; # No, first word goes to left half; Move T,Last; # Yes, load leftover; Movem T,@A; # put it in output; Jrst RH; # and start loading at right half; LH: Caml D,De; # Data emptied?; Jrst LDone; # Yes, return; Move T,@D; # Get current output word; LSH T,16; # Left justify it in 32 rightmost bits; Movem T,@A; # Deposit into Data array; Addi D,1; # and advance output pointer; RH: Caml D,De; # Data emptied?; Jrst RDone; # Yes, return; Move T,@D; # Get current output word; Andi T,'177777; # Mask off possible sign bits; Iorm T,@A; # Put it into output word; Addi D,1; # Advance output pointer; Aoja A,LH; # Advance input pointer; LDone: Setzm LastX; # No leftovers; Jrst DONE; # So bye bye; RDone: Setom LastX; # Have a leftover; Move T,@A; # Get it; Movem Last; # and save for next call (CAN LOSE LAST SAMPLE); Subi A,1; # Back up one; DONE: Sub A,OutArr; # Number of output samples; Movem A,Nout; END "O16c1"; ARRYOUT(Chan,OutArr[1],Nsamps); END "W16B1C"; # Ridiculous conversion routines; # These are necessary when input or output packing is not 16-bit. They should be eliminated when SANDI,SANDO support Xpack. ; PROCEDURE I16toR(INTEGER ARRAY Y; REAL ARRAY Yr; INTEGER N); # Convert 16-bit integer to Real; BEGIN INTEGER i; REAL Scl; Scl _ 2^-15; FOR i_0 Thru N-1 DO Yr[i] _ Y[i]*Scl; RETURN; END; INTEGER PROCEDURE RtoI16(REAL ARRAY Yr; INTEGER ARRAY Y; INTEGER N); # Convert Real to 16-bit integer; BEGIN INTEGER i; INTEGER Bust; REAL Scl; Scl _ 2^15-1; Bust _ 0; FOR i_0 Thru N-1 DO BEGIN REAL Yri; Yri _ Yr[i]; IF ABS(Yri)>1 THEN BEGIN Bust_Bust+1; Yri_Yri/ABS(Yri); END; Y[i] _ Yri*Scl; END; RETURN(Bust); END; # Boolean, integer, and real input routines with help facility; REDEFINE #="comment",thru=" step 1 until ",crlf="('15)&('12)",ALT="'175"; EXTERNAL INTEGER !skip!; BOOLEAN PROCEDURE WantH(BOOLEAN Default; STRING Help,Prompt); BEGIN STRING Ttystr; BOOLEAN Yes,No,Val; LABEL p; p: PRINT(Prompt,"?",(IF Help THEN " (Type ? for help)" ELSE NULL),":"); Ttystr _ INCHRW; IF Ttystr="?" AND Help THEN BEGIN PRINT(CrLf,Help,CrLf); GO TO P; END; IF !skip!=ALT THEN PRINT('15&'12); # ALT means no anywhere; IF Ttystr = CR THEN Ttystr _ INCHRW; # Flush LF in CRLF; Yes _ ((Ttystr="Y") OR (Ttystr="y")); No _ ((Ttystr="N") OR (Ttystr="n")); Val _ (Yes OR (Default AND NOT No)); PRINT((IF VAL THEN " Yes" ELSE " No"),CrLf); RETURN(Val); END; BOOLEAN PROCEDURE RealInH(REFERENCE REAL X; STRING Help,Title); BEGIN "RealInH" BOOLEAN Num; STRING t; INTEGER i; LABEL p; p: PRINT(Title," = ",X," _ ",(IF Help THEN "(Type ? for help) " ELSE NULL)); !skip!_0; t_INCHWL; IF t="?" AND Help THEN BEGIN PRINT(Help,CrLf); GO TO P; END; IF !skip!=ALT THEN BEGIN PRINT(CRLF); RETURN(FALSE); END; Num_FALSE; FOR i_1 Thru LENGTH(t) DO IF "0" LEQ t[i FOR 1] LEQ "9" THEN BEGIN Num_TRUE; DONE; END; IF Num THEN PRINT(Title," set to ",x_REALSCAN(t,0),CRLF) ELSE PRINT(Title, " not changed",CRLF); RETURN(TRUE); END "RealInH"; BOOLEAN PROCEDURE IntInH(REFERENCE INTEGER X; STRING Help,Title); BEGIN "IntInH" BOOLEAN Num; STRING t; INTEGER i; LABEL p; p: PRINT(Title," = ",X," _ ",(IF Help THEN "(Type ? for help) " ELSE NULL)); !skip!_0; t_INCHWL; IF t="?" AND Help THEN BEGIN PRINT(Help,CrLf); GO TO P; END; IF !skip!=ALT THEN BEGIN PRINT(CRLF); RETURN(FALSE); END; Num_FALSE; FOR i_1 Thru LENGTH(t) DO IF "0" LEQ t[i FOR 1] LEQ "9" THEN BEGIN Num_TRUE; DONE; END; IF Num THEN PRINT(Title," set to ",x_INTSCAN(t,0),CRLF) ELSE PRINT(Title, " not changed",CRLF); RETURN(TRUE); END "IntInH"; # Declarations; INTEGER ARRAY OH[0:128]; REAL ARRAY Head[1:128]; INTEGER Ny,Nx,LpScl,Nin,Xoff,Brk,Tlen,Nmult,Nwing,OutPack; REAL Maxamp,Tmp,Factor,Froll,Beta,RtScl; BOOLEAN FAIL,Eof,NewFilter,InterpFilt,FastInput,FastOutput; STRING s,Ifile,Idev,Ofile,Odev; STRING FactorHelp,NmultHelp,FrollHelp,BetaHelp,InterpFiltHelp; LABEL P1,P2,P3,P4; INTERNAL RECORD!CLASS File ( INTEGER Channel, #chans, Pack, Spw, #samps, #elements, #ticks, Flag, State, Type; # 0 = sound file, 1 = fft file ; REAL Clock, Maxamp, # n dB down if fft file ; Bthresh, # Begin of segment threshold ; Ethresh, # End of segment threshold ; SegDur, # Minimum segment duration ; GapDur; # Minimum duration between segments ; STRING Dev, Name, Mode, Packstr, Text ); RECORD!POINTER (File) InpFptr, OutFptr; DEFINE InpF(x) = {File:x[InpFptr]}, OutF(x) = {File:x[OutFptr]}; EXTERNAL PROCEDURE Rhead ( RECORD!POINTER (File) Rp; REFERENCE REAL ARRAY Hdbuf; REFERENCE INTEGER FAIL; BOOLEAN Silence (FALSE) ); # Rhead (rp, hdbuf, Fail, true if don't want printout); EXTERNAL PROCEDURE Whead ( REFERENCE REAL ARRAY Head; REAL Clock; INTEGER Pack, #chans, #ticks; REAL Maxamp; STRING Cstr ); # Whead (head, clock, pack, #chans, #ticks, maxamp, cstr); # Startup; PRINT(CrLf,"Sampling rate conversion program ", COMPILER!BANNER[LENGTH(SCANC(COMPILER!BANNER,Tab,"","sinz"))+11 FOR 17],CrLf); Trpini('14); COMMENT Enable overflow interrupts; SETFORMAT(0,6); InpFptr _ NEW!RECORD (File); OutFptr _ NEW!RECORD (File); OutPack _ Sam16; # Samson box 16-bit format; P1:DO BEGIN OUTSTR("For mono files (R EXTRAC to split out a channel)."&Crlf& "Input "&Def_snd_in_ext&" file: "); Ifile_Default(INCHWL,Idev,Def_snd_in_ext,Def_snd_in_dev); OUTSTR(Tab&"Reading file "&Idev&":"&Ifile&Crlf); OPEN(InpF(Channel)_GETCHAN,Idev,'17,0,0,200,Brk,Eof); LOOKUP(InpF(Channel),Ifile,FAIL); IF FAIL THEN PRINT("File not found",CrLf); END UNTIL NOT FAIL; Rhead(InpFptr,Head,FAIL,TRACE); # Read and print header; Tlen_InpF(#samps); OUTSTR(Crlf&"Output "&Def_snd_out_ext&" file: "); Ofile_Default(INCHWL,Odev,Def_snd_out_ext,Def_snd_out_dev); OUTSTR(Tab&"Writing file "&Odev&":"&Ofile&Crlf); OPEN(OutF(Channel)_GETCHAN,Odev,'17,0,0,200,Brk,Eof); ENTER(OutF(Channel),Ofile,FAIL); IF FAIL THEN USERERR(0,0,"Can't ENTER output file"); FactorHelp _ " Factor is the amount by which the sampling rate is changed. If the sampling rate of the input signal is Srate1, then the sampling rate of the output is Factor*Srate1. "; Factor _ 44100/InpF(Clock); # Default to SONY F1 rate; P2: IF NOT RealInH(Factor,FactorHelp,"Sampling-rate conversion factor") THEN GO TO P1; Maxamp_0; s _ "From "&Idev&":"&Ifile&Crlf& "Comment = "&InpF(Text)&Crlf& "Sampling rate conversion by factor of "&CVF(Factor)&Crlf& "NOT FINISHED YET"&Crlf&Crlf; OH[0]_0; OutF(Dev ) _ Odev; OutF(Name ) _ Ofile; OutF(Clock ) _ InpF(Clock)*Factor; OutF(#chans) _ InpF(#chans); # better be 1 or 2; OutF(#ticks) _ InpF(#ticks)/Factor+0.5; # Unnecessary field - grr; OutF(#samps) _ 0; OutF(Text ) _ s; P3: IF NOT IntInH(OutF(Pack)_OutPack," Packing codes are (0=12b, 1=18b, 3=FP, 4=16b SAM, 5=20b SAM) ","Output packing code") THEN GO TO P2; # Write current version of header (in head) in output file. ; Whead(Head,OutF(Clock),OutF(Pack),OutF(#chans),OutF(#ticks),0,OutF(Text)); USETO(OutF(Channel),1); ARRYOUT(OutF(Channel),Head[1],128); # IO buffer sizes; Nx_TwoTracks; # Two disk tracks; Ny_Nx*Factor+0.5+1; # Round and add 1 for safety; CALL (InpF(Channel), "SHOWIT"); # Simulate an2x for the input channel; # FastInput _ (InpF(Pack)=Sam16); # FastOutput _ (OutF(Pack)=Sam16); # IF NOT FastInput THEN PRINT(" This program runs faster on 16-bit data",CrLf); # Thanks to SANDIO; FastInput _ FastOutput _ FALSE; # Buggy; # Conversion constants; # Old: REQUIRE "SRC.PRM" SOURCE!FILE; # Set converter parameters; # Conversion constants read by SRC.SAI[SYS,MUS], INTRP.SAI[XF,JOS] et al; # ************************************************************ Do not make the macro constants on this page into variables because they are needed at compile time to generate the QUICK!CODE in sampling rate conversion subroutines. ************************************************************ ; DEFINE Nhc = 9, Na = 8, Nh = 16, Nb = 16, Nhxn = 14, NLpScl = 17, Nhg = Nh-Nhxn, Np = Nhc+Na, Npc = 2^Nhc, Pmask = 2^Np-1, Amask = 2^Na-1; # Npc is the number of look-up values available for the lowpass filter between the beginning of its impulse response and the "cutoff time" of the filter. The cutoff time is defined as the reciprocal of the lowpass-filter cutoff frequency in Hz. For example, if the lowpass filter were a sinc function, Npc would be the index of the impulse-response lookup-table corresponding to the first zero-crossing of the sinc function. (The inverse first zero-crossing time of a sinc function equals its nominal cutoff frequency in Hz.) Npc must be a power of 2 due to the details of the current implementation. The default value of 512 is sufficiently high that using linear interpolation to fill in between the table entries gives approximately 16-bit accuracy in filter coefficients. Nhc is log base 2 of Npc. Na is the number of bits devoted to linear interpolation of the filter coefficients. Np is Na + Nhc, the number of bits to the right of the binary point in the integer "time" variable. To the left of the point, it indexes the input array (X), and to the right, it is interpreted as a number between 0 and 1 sample of the input X. Np must be less than 18 in this implementation. Nh is the number of bits in the filter coefficients. The sum of Nh and the number of bits in the input data (typically 16) cannot exceed 36. Thus Nh should be less than 20 (20 does not work). The largest filter coefficient should nearly fill 16 bits (32767). Nb is the number of bits in the input data. The sum of Nb and Nh cannot exceed 36. Nhxn is the number of bits to right shift after multiplying each input sample times a filter coefficient. It can be as great as Nh and as small as 0. Nhxn = Nh-2 gives 2 guard bits in the multiply-add accumulation. If Nhxn=0, the accumulation will soon overflow 36 bits. Nhg is the number of guard bits in mpy-add accumulation ( = Nh-Nhxn). NLpScl is the number of bits allocated to the unity-gain normalization factor. The output of the lowpass filter is multiplied by LpScl and then right-shifted NLpScl bits. To avoid overflow, we must have Nb+Nhg+NlpScl < 36. ; IFC NOT (Np<18) THENC REQUIRE " *** FATAL *** Np<18 DOES NOT HOLD" MESSAGE; ENDC # IF NOT Np<18 THEN USERERR(0,0, " Page 10 of SRC.SAI[SYS,MUS] has violated Np<18"); IFC NOT (Nb+Nhg+NlpScl<36) THENC REQUIRE " *** FATAL *** Nb+Nhg+NlpScl<36 DOES NOT HOLD" MESSAGE; ENDC # IF NOT Nb+Nhg+NlpScl<36 THEN USERERR(0,0, " Page 10 of SRC.SAI[SYS,MUS] has violated Nb+Nhg+NlpScl<36"); IFC NOT (Nh+Nb<36) THENC REQUIRE " *** FATAL *** Nh+Nb<36 DOES NOT HOLD" MESSAGE; ENDC # IF NOT Nh+Nb<36 THEN USERERR(0,0, " Page 10 of SRC.SAI[SYS,MUS] has violated Nh+Nb<36"); # Set up conversion parameters; NmultHelp _ " Nmult is the length of the symmetric FIR lowpass filter used by the sampling rate converter. It must be odd. This is the number of multiplies per output sample for up-conversions (Factor>1), and is the number of multiplies per input sample for down-conversions (Factor<1). Thus if the rate conversion is Srate2 _ Factor*Srate1, then you have Nmult*Srate1*(Factor MAX 1) multiplies per second of real time. Naturally, higher Nmult gives better lowpass-filtering at the expense of longer compute times. Nmult should be odd because it is the length of a symmetric FIR filter, and the current implementation requires a coefficient at the time origin. "; Nmult _ 13; FrollHelp _ " Froll determines the frequency at which the lowpass filter begins to roll-off. If Froll=1, then there is no ""guard zone"" and the filter roll-off region will be aliased. If Froll=0.85, for example, then the filter begins rolling off at 0.85*Srate/2, so that by Srate/2, the filter is well down and aliasing is reduced. Since aliasing distortion is worse by far than loss of the high-frequency spectral amplitude, Froll<1 is highly recommended. The default of .85 sacrifices the upper 15% of the spectrum as an anti-aliasing guard zone. "; Froll _ 0.85; BetaHelp _ " Beta trades the rejection of the lowpass filter against the transition width from passband to stopband. Larger Beta means a slower transition and greater stopband rejection. See Rabiner and Gold (Th. and App. of DSP) under Kaiser windows for more about Beta. The following table from Rabiner and Gold gives some feel for the effect of Beta: All ripples in dB, width of transition band = D*N where N = window length BETA D PB RIP SB RIP 2.120 1.50 +-0.27 -30 3.384 2.23 .0864 -40 4.538 2.93 .0274 -50 5.658 3.62 .00868 -60 6.764 4.32 .00275 -70 7.865 5.0 .000868 -80 8.960 5.7 .000275 -90 10.056 6.4 .000087 -100 "; Beta _ 5.7; InterpFiltHelp _ " By electing to interpolate the filter look-up table, execution is slowed down but the quality of filtering is higher. Interpolation results in twice as many multiply-adds per sample of output. "; InterpFilt _ TRUE; InterpFilt _ WantH(InterpFilt,InterpFiltHelp,"Interpolate filter"); IF (NewFilter _ WantH(NULL,NULL,"Design new filter")) THEN WHILE TRUE DO BEGIN "ol" REAL t; DO IF NOT IntInH(Nmult,NmultHelp," (Odd) Filter length `Nmult' ") THEN GO TO P3 UNTIL Nmult MOD 2 = 1 AND Nmult>0; IF NOT RealInH(Froll,FrollHelp,"Normalized roll-off frequency (0:1) ") THEN CONTINUE "ol"; IF NOT (0 1; Andi a,Amask; # Interpolation factor between IR samples; Floop: Caml Hp,He; # Check for end of filter; Popj p,; # Done; Move t,@Hp; # Put coeff into t; Jumpe IF,NoFi; # May skip filter interpolation; Move 0,@Hdp; # and load difference into AC0; Imul 0,a; # Multiply by interpolation factor; ASH 0,-Na; # a is logically between 0 and 1; Add t,0; # t is now interpolated filter coefficient; Addi Hdp,Npc; # For next time; NoFi: Imul t,@Xp; # Multiply coeff by the input sample; ASH t,-Nhxn; # Leave some guard bits but come back some; Add v,t; # Filter output; Addi Hp,Npc; # IR step; Add Xp,Inc; # Input signal step. NO CHECK MADE FOR ARRAY BOUNDS; Jrst Floop; DONE: Sub Yp,Y; Movem Yp,1; # Return number of output samples; Move t,Regs; # Restore AC's 2 and up; Addi t,2; Hrlz 0,t; Addi 0,2; Blt 0,'17; END "srcu"; END "SrcUp"; # Sampling rate conversion subroutine; INTEGER PROCEDURE SrcUD(INTEGER ARRAY X,Y; REAL Factor; REFERENCE INTEGER Time; INTEGER Nx; INTEGER Nwing,LpScl; INTEGER ARRAY Imp,ImpD; BOOLEAN Interp(TRUE)); BEGIN "SrcUD" REAL Dh; # Step through filter impulse response; REAL Dt; # Step through input signal; INTEGER EndTime; # When Time reaches EndTime, return to user; INTEGER Dhb,Dtb; # Fixed-point versions of Dh,Dt; Dt _ 1/Factor; # Output sampling period; Dtb _ Dt*2^Np+0.5; # Fixed-point representation; Dh _ Factor*Npc MIN Npc; # Filter sampling period; Dhb _ Dh*2^Na+0.5; # Fixed-point representation; EndTime _ Time + Nx*2^Np; # FAIL coded segment; QUICK!CODE "srcf" LABEL Lwing,Rwing,DONE,Floop,Shift,NoFi,Filter; INTEGER XpS; INTEGER ARRAY Regs[0:'17]; DEFINE Xp=1,Yp=2,a=3,Ph=4,v=5,t=6,Xl=7, Hp='10,Hdp='11,He='13,IF='14,Inc='15,Ho='16,P='17; Hrrz 1,Regs; # Save AC's; Move t,Regs; Addi t,'17; Blt 1,(t); Move Yp,Y; # Output is to first sample of output array; Move Xl,X; # Pointer to beginning of input signal; Move He,Nwing; # Max offset + 1 in filter table; Add He,Imp; Move IF,Interp; # Nonzero when interpolating filter table; Lwing: Setzm v; # Clear multiply-add accumulator; Move Xp,Time; # Get current time corresponding to "now" in X; Caml Xp,EndTime;# Check if all Done; Jrst DONE; LSH Xp,-Np; # Pointer to current input sample; Add Xp,Xl; Movem Xp,Xps; # Save for for right-wing of filter (left wing first); Move Ph,Time; # Absolute time from Xl; Andi Ph,Pmask; # Position between input samples; Setom Inc; # X increment is backwards; Pushj P,Filter; # Perform left wing inner product; Rwing: Move Xp,Xps; # Get back to "now" input sample; Addi Xp,1; # Bump forward 1; Setcm Ph,Ph; # This is offset of corr. coeff. in right wing; Andi Ph,Pmask; # Mask off leading 1's; Movei Inc,1; # Now going forward through X; Pushj P,Filter; # Perform right wing inner product; Imul v,LpScl; # Normalize for unity filter gain; ASH v,-NLpScl; Trne v,2^(Nhg-1); # Need to round?; Addi v,2^(Nhg-1); # Yup; ASH v,-Nhg; # Remove guard bits; Movem v,(Yp); # Deposit output; Addi Yp,1; # NO CHECK ON OUTPUT ARRAY BOUNDS; Move t,Dtb; # Time increment; Addm t,Time; # Move to next sample; Jrst Lwing; # Do Next sample; Filter: Move Ho,Ph; # Np bits specifying time since last X sample; Imul Ho,Dhb; # Get (interpolated) offset into filter table; LSH Ho,-Np; # Offset in (interpolated) table; Floop: Move Hp,Ho; # Place offset into Hp; LSH Hp,-Na; # Shift off interpolation bits to get integer part; Add Hp,Imp; # Add in base address of table; Caml Hp,He; # If we have overrun the table; Popj P,; # return happily; Move t,@Hp; # Get IR sample; Jumpe IF,Nofi; # If interpolating filter IR table; Move Hdp,Ho; # Compute integer part of IR address as before; LSH Hdp,-Na; # But this time; Add Hdp,ImpD; # Access the difference table H[n+1]-H[n]; Move 0,@Hdp; # (For linear interpolation); Move a,Ho; # Get interpolation bits; Andi a,Amask; # which are the lower Na bits; Imul 0,a; # Multiply by interpolation factor; ASH 0,-Na; # a is logically between 0 and 1; Add t,0; # t is now interpolated filter coefficient; NoFi: Imul t,@Xp; # Multiply coeff by the input sample; ASH t,-Nhxn; # Leave some guard bits but come back some; Add v,t; # Filter output; Add Ho,Dhb; # IR step; Add Xp,Inc; # Input signal step. NO CHECK MADE FOR ARRAY BOUNDS; Jrst Floop; DONE: Sub Yp,Y; Movem Yp,1; # Return number of output samples; Move t,Regs; # Restore AC's 2 and up; Addi t,2; Hrlz 0,t; Addi 0,2; Blt 0,'17; END "srcf"; END "SrcUD"; # Set up lowpass filter. If NewFilter is FALSE, we attempt to read the interpolation lowpass filter from a standard disk file. Otherwise we design the lowpass filter (a Kaiser-windowed sinc function) here. ; IF NOT NewFilter THEN BEGIN "Try to find the right disk file" INTEGER Fchan,Brk; BOOLEAN Eof,FAIL; STRING Fname; Fname _ "F"&CVS(Nmult)&"T"&CVMS(Nhc)&".SRC[1,JOS]"; # Construct filename; # Original default filename is "F13T9.SRC[1,JOS]"; OUTSTR(Tab&"Reading file DSK:"&Fname&Crlf); OPEN(Fchan_GETCHAN,"DSK",'17,0,0,200,Brk,Eof); LOOKUP(Fchan,Fname,FAIL); IF FAIL THEN BEGIN PRINT("File not found",CrLf); NewFilter _ TRUE; END ELSE BEGIN "ReadFilterFile" # Filter file format: The first word of the file is (scaleFactor,Length), then there are 127 zeros (to finish a record), then there are "Length" 16-bit impulse response values in the right wing of the impulse response. (Length=Npc*(Nmult+1)/2+1 where originally Npc=2^9, and Nmult=13.) At the next record boundary (one record = 128 words, one word = 36 bits [haha]), there is a similar table of impulse-response successive differences: Diff[i]=Imp[i+1]-Imp[i]. ; INTEGER Wd; Wd _ WORDIN(Fchan); # First word of file is [ScaleFactor,,Length]; LpScl _ Wd LSH -18; # (Integer) Scale factor has to be positive; Nwing _ Wd LAND '777777; # Mask off table length; ARRYIN(Fchan,Imp[0],Nwing); # Array of 16-bit filter coefficients; ARRYIN(Fchan,ImpD[0],Nwing);# Array of 16-bit filter coeff. differences; # Account for increased filter gain when using Factors less than 1; IF Factor<1 THEN LpScl _ LpScl*Factor + 0.5; END "ReadFilterFile"; RELEASE(Fchan); END "Try to find the right disk file"; IF NewFilter THEN BEGIN "Make lowpass filter" REAL DCgain,Scl,Maxh; INTEGER i,Dh; PRINT(" Setting up lowpass filter . . . "); # FirLpf is in JAMLIB. It designs a KaiserWindow*Sinc lowpass filter; FirLpf(ImpR,Nwing,.5*Froll,Npc,Beta); # Compute the DC gain of the lowpass filter, and its maximum coefficient magnitude. Scale the coefficients so that the maximum coefficient just fits in Nh-bit fixed-point, and compute LpScl as the NLpScl-bit (signed) scale factor which when multiplied by the output of the lowpass filter gives unity gain. LpScl is stored for the case Factor GEQ 1. ; DCgain_0; Dh _ Npc; # Filter sampling period for factors greater than or equal to 1; FOR i_Dh STEP Dh UNTIL Nwing-1 DO DCgain_DCgain+ImpR[i]; DCgain_2*DCgain+ImpR[0]; # DC gain of real coefficients; FOR i_0 STEP 1 UNTIL Nwing-1 DO Maxh _ Maxh MAX ABS(ImpR[i]); # ImpR[0] typ.; Scl _ (2^(Nh-1)-1)/Maxh; # Map largest coefficient to 16-bit maximum; LpScl _ ABS(2^(NLpScl+Nh)/(DcGain*Scl)); IF LpScl>(2^18-1) THEN USERERR(0,1, " Oops! Filter scale factor does not fit in halfword."& " Reduce NlpScl or try another filter"); # Scale filter coefficients for Nh bits and convert to integer ; IF ImpR[0]<0 THEN Scl_-Scl; # Need positive first value for LpScl storage; FOR i_0 STEP 1 UNTIL Nwing-1 DO ImpR[i] _ ImpR[i]*Scl; FOR i_0 Thru Nwing-1 DO Imp[i] _ ImpR[i] + 0.5; # Round it; # ImpD makes linear interpolation of the filter coefficients faster; FOR i_0 Thru Nwing-2 DO ImpD[i] _ ImpR[i+1] - ImpR[i] + 0.5; ImpD[Nwing-1] _ Imp[Nwing-1]; # Last coeff. not interpolated; PRINT("Done.",CrLf); BEGIN "write" INTEGER PROCEDURE CALLI(INTEGER Type;INTEGER AC1(0)); START!CODE MOVE 1,AC1; CALLI 1, @Type; END; STRING PROCEDURE PPN; RETURN(CVXSTR(CALLI('400071,0))); IF EQU(PPN," 1JOS") THEN BEGIN "Write lowpass filter" INTEGER Fchan,Brk; BOOLEAN Eof,FAIL; STRING Fname; Fname _ "F"&CVS(Nmult)&"T"&CVMS(Nhc)&".SRC[1,JOS]"; OUTSTR(Tab&"Writing file DSK:"&Fname&Crlf); OPEN(Fchan_GETCHAN,"DSK",'17,0,0,200,Brk,Eof); ENTER(Fchan,Fname,FAIL); IF FAIL THEN PRINT("Can't enter file",CrLf) ELSE BEGIN INTEGER Wd; # Imp[0] _ (Imp[0] LAND '777777) LOR (LpScl LSH 18); Wd _ (LpScl LSH 18) + Nwing; WORDOUT(Fchan,Wd); ARRYOUT(Fchan,Imp[0],Nwing); # Imp[0] _ (Imp[0] LAND '777777); ARRYOUT(Fchan,ImpD[0],Nwing); END; RELEASE(Fchan); END "Write lowpass filter"; END "write"; CkFIR(Imp,Nwing,Factor*Npc MIN Npc,InterpFilt); END "Make lowpass filter" ELSE IF TRACE THEN CkFIR(Imp,Nwing,Factor*Npc MIN Npc,InterpFilt); IF TRACE THEN PRINT(" LpScl = ",LpScl,CrLf); # Inner Loop; Xoff _ ((Nmult+1)/2)*(1/Factor MAX 1); # Reach of lowpass filter wing; Xoff _ Xoff + 10; # Give it some creeping room; Nin _ Nx + 2*Xoff; # Number of points to read each buffer; SETFORMAT(0,2); # For buffer begin-time print-out; IF InpF(#chans)=1 THEN BEGIN "mono" INTEGER Yp,Xp,Time,Ncreep; Saf INTEGER ARRAY X[0:Nin-1],Y[0:Ny-1]; Saf REAL ARRAY Xr[0:Nin-1],Yr[0:Ny-1]; Yp _ 0; # Current sample and length of output file; Xp _ 0; # Current "now"-sample pointer for input; Time _ Xoff LSH Np; # Current-time pointer for converter; RtScl _ .001; # Convert milliseconds to seconds; PRINT("At time (Itime, Xtime, Otime): "); DO BEGIN "RI" INTEGER RunTime,Nout,Novfl; PRINT(CVF(Xp/InpF(Clock))," "); IF Xp-Xoff+Nin>Tlen THEN BEGIN Nx _ Nx - (Xp-Xoff+Nin-Tlen); Nin _ Nx + 2*Xoff; END; RunTime _ CALL(0,"RUNTIM"); IF FastInput THEN I16B1C(InpF(Channel),Xp-Xoff,Nin,X) ELSE BEGIN Sandi(InpF(Channel),Xp-Xoff,Nin,Xr,Tlen,TRUE,InpF(Pack)); Novfl _ RtoI16(Xr,X,Nin); IF Novfl THEN PRINT(" Overflowed 16 bits on input ",Novfl," times.",CrLf); END; PRINT("(",CVF((RunTime _ CALL(0,"RUNTIM") - RunTime)*RtScl),","); IF Factor GEQ 1 THEN Nout _ SrcUp(X,Y,Factor,Time,Nx,Nwing,LpScl,Imp,ImpD,InterpFilt) ELSE Nout _ SrcUD(X,Y,Factor,Time,Nx,Nwing,LpScl,Imp,ImpD,InterpFilt); PRINT(CVF((RunTime _ CALL(0,"RUNTIM") - RunTime)*RtScl),","); IF Trace THEN BEGIN PRINT(CrLf," Time='",CVOS(Time)," = ",Time LSH -Np," LSH Np"); PRINT(" Nout =",Nout); END; Time _ Time - (Nx LSH Np); # Move converter Nx samples back in time; Xp _ Xp + Nx; # Advance by number of samples processed; Ncreep _ (Time LSH -Np) - Xoff; IF Trace THEN PRINT(" Crept ",Ncreep,CrLf); IF Ncreep NEQ 0 THEN BEGIN Time _ Time - (Ncreep LSH Np); # Remove time accumulation; Xp _ Xp + Ncreep; # and add it to read pointer; END; IF Nout>Ny THEN USERERR(0,1,"Output array overflow"); RunTime _ CALL(0,"RUNTIM"); IF FastOutput THEN W16B1C(OutF(Channel),Nout,Y) ELSE BEGIN I16toR(Y,Yr,Nout); IF OutF(Pack)=4 THEN BEGIN INTEGER i; FOR i_0 Thru Nout-1 DO IF ABS(Yr[i])>1 THEN BEGIN Novfl_Novfl+1; Yr[i]_Yr[i]/ABS(Yr[i]); END; IF Novfl>0 THEN PRINT(" Output Overflows 16 bits. Clipped ",Novfl," times.",CrLf); END; Sando(OutF(Channel),Yp,Nout,Yr,OH,Yp,TRUE,OutF(Pack)); END; PRINT(CVF((CALL(0,"RUNTIM") - RunTime)*RtScl),") "); END "RI" UNTIL Xp GEQ Tlen-Xoff; IF NOT FastOutput THEN Sando(OutF(Channel),0,0,Yr,OH,Yp,TRUE,OutF(Pack)); END "mono" ELSE BEGIN "oops" PRINT(Crlf&"Sorry, can only do 1-channel files."); CALL(0,"EXIT"); END "oops"; END "ALAR"; # update header; s_"From "&Idev&":"&Ifile&Crlf& "Comment = "&InpF(Text)&Crlf& "Sampling rate conversion by factor of "&CVF(Factor)&Crlf; Whead (Head, OutF(Clock), OutF(Pack), OutF(#chans), OutF(#ticks), Maxamp,s); USETO (OutF(Channel), 1); ARRYOUT (OutF(Channel), Head [1], 128); RELEASE(OutF(Channel)); END "SRC".