JOS Home
Original source code for resample
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 an 2x 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 (01;
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".