*----------------------------------------------------------------------- * spectr.f - read digital audio, draw spectrogram as a Portable Grey Map. * by Andy Allinger, 2016-2021, released to the public domain * * Permission to use, copy, modify, and distribute this software and * its documentation for any purpose and without fee is hereby * granted, without any conditions or restrictions. This software is * provided "as is" without express or implied warranty. * * OPTIONS * --tune x transpose input up by x semitones * --swab reverse byte order of the input *----------------------------------------------------------------------- PROGRAM SPECTR IMPLICIT NONE * Constants INTEGER BAND, DUR, FLIM, HIRATE, HLIM, LLIM, M, MAXBRT, MNNSUB, $ NLIM, PXPS, VERT REAL BIG, FMARK, FMNNLO, FMNNHI, FSPAC, HPASS, HSTOP, LIT, $ OLAP, PI, PI2, PIP, RESPAR, SIG, TARGN, TOLT *----------------------------------------------------------------------- * Not all of these parameters are configurable. *----------------------------------------------------------------------- PARAMETER (BAND = 32, ! number of transform bands $ BIG = 1.E+36, ! arbitrary large number $ DUR = 5 * 60, ! max file duration [seconds] $ FLIM = 128, ! most factors $ HIRATE = 48000, ! fast sample rate $ FMNNLO = 52.5, ! lowest frequency as MIDI note number $ FMNNHI = 84.5, ! highest frequency $ LIT = 1.E-36, ! very small number $ MAXBRT = 255, ! maximum value in a .pgm file $ M = 17, ! multiple of fundemental frequency $ MNNSUB = 52) ! MIDI note below transform range PARAMETER (NLIM = 27720, ! largest window size $ OLAP = 6., ! 1/overlap of successive windows $ PI = 3.141592654, $ PI2 = 9.869604401, ! pi^2 $ PIP = 0.035, ! a brief duration $ PXPS = 240, ! Pixels per second $ RESPAR = 1. / 1024., ! desired precision: -60 dB $ SIG = 2., ! Lanczos parameter $ TARGN = 3000., ! target value for N $ TOLT = 17.3123405, ! twelve over log two $ VERT = 640) ! height of output image PARAMETER ( $ FSPAC = (FMNNHI - FMNNLO) / VERT, ! spacing between pixels $ HPASS = 1. - RESPAR, ! passband transmission $ HSTOP = RESPAR, ! stopband transmission $ HLIM = DUR * PXPS, ! max image width $ LLIM = DUR * HIRATE) ! max # samples PARAMETER (FMARK = FMNNLO - 0.5 * FSPAC) ! center freq of 0th pixel * Variables INTEGER $ I, J, K, ! counters $ FACTR(FLIM), ! factors of RATE $ HORZ, ! image width $ IERR, ! error code $ ILO, IHI, ! coefficient indices $ IND(BAND), ! permutation vector $ L, ! length of input, samples $ N, ! window size $ NEF(HLIM), ! number of transforms each pixel $ R, RATE, ! sample rate, input sample rate $ STRIDE, ! factor to downsample by $ TARGR ! target sampling rate REAL $ A, B, ! bounds of interpolation $ AMPL(HLIM,VERT), ! intensity per frequency, time $ DIF, EMIN, EMAX, ! a difference, extreme values $ F, ! frequency $ FPASS, FSTOP, ! passband, stopband frequencies $ FSAMP, ! sample rate as real $ HOP, ! # samples to advance $ RAXBRT, ! MAXBRT as REAL $ S(LLIM), SF(LLIM), ! samples, filtered samples $ TOT, ! a total $ TRNPOS, ! amount to transpose $ U, ! a scale factor $ FREQ(HLIM,BAND), ! frequency of the ith coefficient $ X(BAND), Y(BAND), ! workspace arrays $ XI ! frequency of pixel DOUBLE PRECISION $ T(2,NLIM), ! transform matrix $ W(2,NLIM), ! window function, its derivative $ Z(NLIM) ! samples in mono LOGICAL BYSW, LEXIST CHARACTER*256 IFIL, OFIL, OPT ! filenames * External functions INTEGER IARGC, LTRIM LOGICAL SAFDIV DATA FACTR / $ 1, 2, 3, 4, 5, 6, 7, 9, $ 10, 12, 14, 15, 18, 20, 21, 25, $ 28, 30, 35, 36, 42, 45, 49, 50, $ 60, 63, 70, 75, 84, 90, 98, 100, $ 105, 126, 140, 147, 150, 175, 180, 196, $ 210, 225, 245, 252, 294, 300, 315, 350, $ 420, 441, 450, 490, 525, 588, 630, 700, $ 735, 882, 900, 980, 1050, 1225, 1260, 1470, $ 1575, 1764, 2100, 2205, 2450, 2940, 3150, 3675, $ 4410, 4900, 6300, 7350, 8820, 11025, 14700, 22050, $ 48 * 44100 / *----------------------------------------------------------------------- * Begin. *----------------------------------------------------------------------- IERR = 0 N = IARGC () IF (N < 2) THEN PRINT *, 'spectr: files not specified' STOP ELSE IF (N > 5) THEN PRINT *, 'spectr: too many arguments!' STOP END IF CALL GETARG (N-1, IFIL) ! read file name from command line CALL GETARG (N, OFIL) * Options processing loop TRNPOS = 0. BYSW = .FALSE. DO I = 1, N-2 CALL GETARG (I, OPT) CALL UCASE (OPT) IF (INDEX(OPT, '--SWAB') > 0) BYSW = .TRUE. IF (INDEX(OPT, '--TUNE') > 0) THEN CALL GETARG(I+1, OPT) READ (UNIT=OPT, FMT=*, IOSTAT=IERR) TRNPOS IF (IERR .NE. 0) PRINT *, 'trouble reading retuning option!' IF (TRNPOS < MNNSUB + BAND - 127 .OR. TRNPOS > MNNSUB) THEN PRINT *, 'transpose out of bounds' STOP END IF END IF END DO *----------------------------------------------------------------------- * Check files for existence *----------------------------------------------------------------------- INQUIRE (FILE=IFIL, EXIST=LEXIST, IOSTAT=IERR) IF (IERR .NE. 0) STOP 'trouble inquiring of input file!' IF (.NOT. LEXIST) THEN L = LTRIM(IFIL) PRINT *, 'File ', IFIL(1:L), ' not found.' STOP END IF INQUIRE (FILE=OFIL, EXIST=LEXIST, IOSTAT=IERR) IF (IERR .NE. 0) STOP 'trouble inquiring of output file!' IF (LEXIST) THEN L = LTRIM(OFIL) PRINT *, 'File: ', OFIL(1:L), ' already exists. ', $ 'Please choose another name.' STOP END IF CALL LFIL (IFIL, N, IERR) IF (IERR .NE. 0) STOP 'trouble getting length of input file' *----------------------------------------------------------------------- * Read in samples *----------------------------------------------------------------------- CALL RWAVMO (IFIL, N, BYSW, S, LLIM, L, RATE, IERR) IF (IERR .NE. 0) THEN PRINT *, 'trouble reading from file: ', IFIL STOP END IF FSAMP = FLOAT(RATE) HORZ = INT (DBLE(L) * DBLE(PXPS) / DBLE(FSAMP)) + 1 ! width of image IF (L < PIP * FSAMP) STOP 'too short a file!' IF (L > LLIM .OR. HORZ > HLIM) STOP 'too big a file!' IF (RATE .EQ. 44100) THEN K = 81 ELSE K = 0 DO J = 1, RATE IF (MOD(RATE, J) .EQ. 0) THEN IF (K .GE. FLIM) STOP 'too many factors of sampling rate!' K = K + 1 FACTR(K) = J END IF END DO END IF *----------------------------------------------------------------------- * Short-time Fourier transforms. M @ R = N @ F *----------------------------------------------------------------------- DO J = 1, BAND F = 55. * 2.**((MNNSUB + J - 33 - TRNPOS) / 12.) TARGR = NINT(TARGN * F / M) CALL IBSGE (FACTR, K, TARGR, I) R = FACTR(I) * Filter. DO I = 1, L SF(I) = S(I) END DO IF (R < RATE) THEN ! lo pass FPASS = F * ((2. * M + 1.) / (2. * M)) FSTOP = 0.5 * R CALL FILTER (SF, L, FSAMP, FPASS, FSTOP, HPASS, HSTOP, IERR) IF (IERR .NE. 0) PRINT *, 'lo pass filter returns #', IERR END IF FPASS = F * ((2. * M - 1.) / (2. * M)) ! hi pass FSTOP = F / M CALL FILTER (SF, L, FSAMP, FPASS, FSTOP, HPASS, HSTOP, IERR) IF (IERR .NE. 0) PRINT *, 'hi pass filter returns #', IERR N = NINT(M * R / F) IF (N > NLIM) THEN PRINT *, 'Window size too big! Try a lower sample rate.' STOP END IF STRIDE = RATE / R HOP = MIN(N * STRIDE / OLAP, FSAMP / PXPS) CALL FRAT (SF, L, M, N, R, STRIDE, HOP, HORZ, $ FREQ(1,J), AMPL(1,J), NEF, T, W, Z, IERR) IF (IERR .NE. 0) PRINT *, 'frat returns #', IERR END DO *----------------------------------------------------------------------- * Convert to MIDI note number. *----------------------------------------------------------------------- DO K = 1, BAND DO J = 1, HORZ FREQ(J,K) = TOLT * LOG(FREQ(J,K) / 55.) + 33. END DO END DO *----------------------------------------------------------------------- * Interpolate. *----------------------------------------------------------------------- DO K = 1, HORZ DO J = 1, BAND ! transpose copy X(J) = FREQ(K,J) Y(J) = AMPL(K,J) END DO CALL ISORTR (X, IND, BAND) CALL RORDER (X, IND, BAND) CALL RORDER (Y, IND, BAND) DO J = 1, VERT ! interpolate XI = FMARK - TRNPOS + FSPAC * J A = XI - FSPAC * SIG B = XI + FSPAC * SIG CALL RBSGT (X, BAND, A, ILO) CALL RBSLT (X, BAND, B, IHI) TOT = 0. DO I = ILO, IHI DIF = (X(I) - XI) / FSPAC A = SIG * SIN(PI * DIF) * SIN(PI * DIF / SIG) ! Lanczos kernel B = PI2 * DIF**2 IF (SAFDIV(A, B)) THEN U = A / B ELSE U = 1. END IF TOT = TOT + Y(I) * U END DO AMPL(K,J) = TOT END DO END DO ! next time step *----------------------------------------------------------------------- * Scale image intensities 0...MAXBRT *----------------------------------------------------------------------- EMIN = +BIG ! find extreme values EMAX = -BIG DO K = 1, VERT DO J = 1, HORZ IF (AMPL(J,K) < EMIN) EMIN = AMPL(J,K) IF (AMPL(J,K) > EMAX) EMAX = AMPL(J,K) END DO END DO * Logarithms U = EMAX * RESPAR IF (U > LIT) THEN DO K = 1, VERT DO J = 1, HORZ AMPL(J,K) = LOG(MAX(0., AMPL(J,K)) + U) END DO END DO EMIN = LOG(MAX(0., EMIN) + U) EMAX = LOG(MAX(0., EMAX) + U) END IF RAXBRT = FLOAT(MAXBRT) DIF = EMAX - EMIN IF (SAFDIV(RAXBRT, DIF)) THEN U = RAXBRT / DIF ELSE U = 1. END IF DO K = 1, VERT DO J = 1, HORZ AMPL(J,K) = (AMPL(J,K) - EMIN) * U AMPL(J,K) = MIN(MAX(0., AMPL(J,K)), RAXBRT) END DO END DO *----------------------------------------------------------------------- * write out *----------------------------------------------------------------------- CALL WPGM (OFIL, AMPL, HLIM, HORZ, VERT, IERR) IF (IERR .NE. 0) PRINT *, 'wpgm returns #', IERR STOP 'program complete' END ! of spectr *----------------------------------------------------------------------- * rwavmo - read a wave file into mono * *___Name_____Type____________In/Out____Description______________________ * IFIL Character*(*) In Name of input file * FILLEN Integer In Length of input file in bytes * BYSW Logical In Reverse byte order of floats * S(KS) Real Out Buffer of samples * KS Integer In Array bound of S * L Integer Out Length of S * RATE Integer Out Samples / second * IERR Integer Out Error code *----------------------------------------------------------------------- SUBROUTINE RWAVMO (IFIL, FILLEN, BYSW, S, KS, L, RATE, IERR) IMPLICIT NONE INTEGER FILLEN, KS, L, RATE, IERR REAL S(KS) LOGICAL BYSW CHARACTER*(*) IFIL * Local variables INTEGER*1 B INTEGER $ BLKALN, ! channels @ bits/second / 8 $ BPS, ! bits/sampl $ BYTRAT, ! channels @ bits/second / 8 $ CHANEL, ! two for stereo $ CHNKSZ, ! samples @ channels @ bits/second / 8 $ I, J, K, ! count bytes, count samples $ IBEG, IEND, ! byte positions $ IOU, ! i/o unit $ ISH, $ NBY, ! bytes per integer $ POS, ! offset in bytes $ SAMPL, $ SHIFT, ! offset for signed integers $ STRIDE, ! count up or down $ THRESH ! max positive value of signed integer REAL FLOSAM, DIVID, FCH LOGICAL FEXIST, FP CHARACTER*4 $ IDRIFF, $ IDRIFX, $ IDWAVE, $ IDFMT, $ IDDATA, $ TAG * External function INTEGER IOUNIT * Constants DATA IDRIFF / 'RIFF' / DATA IDRIFX / 'RIFX' / DATA IDWAVE / 'WAVE' / DATA IDFMT / 'fmt ' / DATA IDDATA / 'data' / *----------------------------------------------------------------------- * Open input file *----------------------------------------------------------------------- IERR = 0 L = 0 RATE = 0 INQUIRE (FILE=IFIL, EXIST=FEXIST, IOSTAT=IERR) IF (IERR .NE. 0) RETURN IF (.NOT. FEXIST) THEN PRINT *, 'file not found: ', IFIL IERR = 1 RETURN END IF IOU = IOUNIT () OPEN (UNIT=IOU, FILE=IFIL, STATUS='OLD', ACCESS='DIRECT', $ RECL=1, FORM='UNFORMATTED', IOSTAT=IERR) IF (IERR .NE. 0) RETURN *----------------------------------------------------------------------- * Read RIFF *----------------------------------------------------------------------- POS = 0 TAG = ' ' DO J = 0, 3 POS = POS + 1 IF (POS > FILLEN) GO TO 90 READ (UNIT=IOU, REC=POS, IOSTAT=IERR) B IF (IERR .NE. 0) THEN PRINT *, 'trouble reading byte ', J, ' of RIFF' RETURN END IF TAG(J+1:J+1) = CHAR(B) END DO IF (TAG .EQ. IDRIFF) THEN STRIDE = +1 ELSE IF (TAG .EQ. IDRIFX) THEN STRIDE = -1 ELSE IERR = 2 RETURN END IF *----------------------------------------------------------------------- * Read file size. *----------------------------------------------------------------------- I = 0 DO J = 0, 3 POS = POS + 1 IF (POS > FILLEN) GO TO 90 READ (UNIT=IOU, REC=POS, IOSTAT=IERR) B IF (IERR .NE. 0) THEN PRINT *, 'trouble reading file size byte ', J RETURN END IF CALL BYTGET (B, I, J) END DO IF (I .NE. FILLEN - 8) THEN PRINT *, 'i is ', I, ' but measured file length is: ', FILLEN IERR = 3 RETURN END IF *----------------------------------------------------------------------- * Read WAVE *----------------------------------------------------------------------- TAG = ' ' DO J = 0, 3 POS = POS + 1 IF (POS > FILLEN) GO TO 90 READ (UNIT=IOU, REC=POS, IOSTAT=IERR) B IF (IERR .NE. 0) THEN PRINT *, 'trouble reading byte ', J, ' of WAVE' RETURN END IF TAG(J+1:J+1) = CHAR(B) END DO IF (TAG .NE. IDWAVE) THEN IERR = 4 RETURN END IF *----------------------------------------------------------------------- * Read chunk descriptor *----------------------------------------------------------------------- 10 TAG = ' ' DO J = 0, 3 POS = POS + 1 IF (POS > FILLEN) GO TO 90 READ (UNIT=IOU, REC=POS, IOSTAT=IERR) B IF (IERR .NE. 0) THEN PRINT *, 'trouble reading byte ', J, ' of fmt' RETURN END IF TAG(J+1:J+1) = CHAR(B) END DO *----------------------------------------------------------------------- * Read 16 *----------------------------------------------------------------------- I = 0 DO J = 0, 3 POS = POS + 1 IF (POS > FILLEN) GO TO 90 READ (UNIT=IOU, REC=POS, IOSTAT=IERR) B IF (IERR .NE. 0) THEN PRINT *, 'trouble reading 16 byte ', J RETURN END IF CALL BYTGET (B, I, J) END DO IF (TAG .NE. IDFMT) THEN POS = POS + I GO TO 10 ELSE IF (I .NE. 16) THEN IERR = 5 RETURN END IF *----------------------------------------------------------------------- * Read 1 for PCM *----------------------------------------------------------------------- ISH = 0 DO J = 0, 1 POS = POS + 1 IF (POS > FILLEN) GO TO 90 READ (UNIT=IOU, REC=POS, IOSTAT=IERR) B IF (IERR .NE. 0) THEN PRINT *, 'Trouble reading PCM byte ', J RETURN END IF CALL BYTGET (B, ISH, J) END DO IF (ISH .EQ. 1) THEN FP = .FALSE. ELSE IF (ISH .EQ. 3) THEN FP = .TRUE. ELSE PRINT *, 'unsupported codec #', ISH, ' (must be PCM)' IERR = 6 RETURN END IF *----------------------------------------------------------------------- * Read # channels *----------------------------------------------------------------------- CHANEL = 0 DO J = 0, 1 POS = POS + 1 IF (POS > FILLEN) GO TO 90 READ (UNIT=IOU, REC=POS, IOSTAT=IERR) B IF (IERR .NE. 0) THEN PRINT *, 'trouble reading # channels byte ', J RETURN END IF CALL BYTGET (B, CHANEL, J) END DO IF (CHANEL < 1 .OR. CHANEL > 2) THEN PRINT *, 'bad # channels!' IERR = 7 RETURN END IF *----------------------------------------------------------------------- * Read sample rate *----------------------------------------------------------------------- RATE = 0 DO J = 0, 3 POS = POS + 1 IF (POS > FILLEN) GO TO 90 READ (UNIT=IOU, REC=POS, IOSTAT=IERR) B IF (IERR .NE. 0) THEN PRINT *, 'trouble reading sample rate byte ', J RETURN END IF CALL BYTGET (B, RATE, J) END DO IF (RATE < 2205 .OR. RATE > 1 000 000) THEN PRINT *, 'sample rate out of bounds!' IERR = 8 RETURN END IF *----------------------------------------------------------------------- * Read byte rate *----------------------------------------------------------------------- BYTRAT = 0 DO J = 0, 3 POS = POS + 1 IF (POS > FILLEN) GO TO 90 READ (UNIT=IOU, REC=POS, IOSTAT=IERR) B IF (IERR .NE. 0) THEN PRINT *, 'trouble reading byte rate, byte ', J RETURN END IF CALL BYTGET (B, BYTRAT, J) END DO *----------------------------------------------------------------------- * Read block align *----------------------------------------------------------------------- BLKALN = 0 DO J = 0, 1 POS = POS + 1 IF (POS > FILLEN) GO TO 90 READ (UNIT=IOU, REC=POS, IOSTAT=IERR) B IF (IERR .NE. 0) THEN PRINT *, 'trouble reading block align byte ', J RETURN END IF CALL BYTGET (B, BLKALN, J) END DO *----------------------------------------------------------------------- * Read bits per sample *----------------------------------------------------------------------- BPS = 0 DO J = 0, 1 POS = POS + 1 IF (POS > FILLEN) GO TO 90 READ (UNIT=IOU, REC=POS, IOSTAT=IERR) B IF (IERR .NE. 0) THEN PRINT *, 'trouble reading bits per sample byte ', J RETURN END IF CALL BYTGET (B, BPS, J) END DO IF (BPS < 1 .OR. MOD(BPS, 8) .NE. 0) THEN PRINT *, 'bad bits per sample: ', BPS IERR = 9 RETURN END IF NBY = BPS / 8 IF (BYTRAT .NE. RATE * CHANEL * NBY) THEN PRINT *, '!!! bytrat: ', BYTRAT, ' rate: ', RATE, $ ' chanel: ', CHANEL, ' bits per second: ', BPS IERR = 10 RETURN END IF IF (BLKALN .NE. CHANEL * NBY) THEN PRINT *, '!! blkaln: ', BLKALN, ' chanel: ', CHANEL, $ ' bits per second: ', BPS IERR = 11 RETURN END IF *----------------------------------------------------------------------- * Read descriptor *----------------------------------------------------------------------- 20 TAG = ' ' DO J = 0, 3 POS = POS + 1 IF (POS > FILLEN) GO TO 90 READ (UNIT=IOU, REC=POS, IOSTAT=IERR) B IF (IERR .NE. 0) THEN PRINT *, 'trouble reading data chunk header byte ', J RETURN END IF TAG(J+1:J+1) = CHAR(B) END DO *----------------------------------------------------------------------- * Read size of data chunk *----------------------------------------------------------------------- CHNKSZ = 0 DO J = 0, 3 POS = POS + 1 IF (POS > FILLEN) GO TO 90 READ (UNIT=IOU, REC=POS, IOSTAT=IERR) B IF (IERR .NE. 0) THEN PRINT *, 'trouble reading size of data chunk byte ', J RETURN END IF CALL BYTGET (B, CHNKSZ, J) END DO IF (TAG .NE. IDDATA) THEN POS = POS + CHNKSZ GO TO 20 END IF *----------------------------------------------------------------------- * Read samples *----------------------------------------------------------------------- L = CHNKSZ / CHANEL / NBY IF (L > KS) THEN PRINT *, 'too big a file!' IERR = 12 RETURN END IF IF (STRIDE > 0) THEN IBEG = 0 IEND = NBY - 1 ELSE IBEG = NBY - 1 IEND = 0 END IF IF (CHANEL .EQ. 1) THEN FCH = 1. ELSE IF (CHANEL .EQ. 2) THEN FCH = 0.5 ELSE PRINT *, 'number of channels: ', CHANEL IERR = 13 RETURN END IF IF (FP) GO TO 30 IF (NBY .EQ. 1) THEN ! unsigned THRESH = 0 SHIFT = 128 DIVID = 1. / 128. ELSE IF (NBY .EQ. 2) THEN ! signed THRESH = 32767 SHIFT = 65536 DIVID = 1. / 32768. ELSE IF (NBY .EQ. 3) THEN ! signed THRESH = 8388607 SHIFT = 16777216 DIVID = 1. / 8388608. ELSE PRINT *, '! nby is ', NBY IERR = 14 RETURN END IF *----------------------------------------------------------------------- * Main loop. *----------------------------------------------------------------------- DO K = 1, L S(K) = 0. DO J = 1, CHANEL SAMPL = 0 DO I = IBEG, IEND, STRIDE POS = POS + 1 IF (POS > FILLEN) GO TO 90 READ (UNIT=IOU, REC=POS, IOSTAT=IERR) B IF (IERR .NE. 0) THEN PRINT *, 'trouble reading byte ', POS, ' of ', IFIL RETURN END IF CALL BYTGET (B, SAMPL, I) ! copy bits END DO IF (SAMPL > THRESH) SAMPL = SAMPL - SHIFT S(K) = S(K) + FLOAT(SAMPL) * DIVID END DO S(K) = S(K) * FCH END DO ! loop on samples GO TO 80 *----------------------------------------------------------------------- * Floating-point read *----------------------------------------------------------------------- 30 CLOSE (UNIT=IOU, IOSTAT=IERR) IF (IERR .NE. 0) THEN PRINT *, 'trouble closing file: ', IFIL RETURN END IF IF (NBY .NE. 4) THEN PRINT *, 'floating point numbers must be 32 bits' IERR = 15 RETURN END IF POS = POS / 4 IOU = IOUNIT () OPEN (UNIT=IOU, FILE=IFIL, STATUS='OLD', ACCESS='DIRECT', $ RECL=4, FORM='UNFORMATTED', IOSTAT=IERR) DO K = 1, L S(K) = 0. DO J = 1, CHANEL POS = POS + 1 IF (POS > FILLEN) GO TO 90 READ (UNIT=IOU, REC=POS, IOSTAT=IERR) FLOSAM IF (IERR .NE. 0) THEN PRINT *, 'trouble reading byte ', POS, ' of ', IFIL RETURN END IF IF (BYSW) CALL SWAB32 (FLOSAM) S(K) = S(K) + FLOSAM END DO S(K) = S(K) * FCH END DO ! loop on samples *----------------------------------------------------------------------- * Done with input file *----------------------------------------------------------------------- 80 CLOSE (UNIT=IOU, IOSTAT=IERR) IF (IERR .NE. 0) THEN PRINT *, 'trouble closing file: ', IFIL RETURN END IF RETURN ! success 90 PRINT *, 'try to read past end of file!' IERR = 16 RETURN END ! of rwavmo *----------------------------------------------------------------------- SUBROUTINE SWAB32 (R) IMPLICIT NONE REAL R INTEGER*1 B0, B1, B2, B3 INTEGER IV REAL RV EQUIVALENCE (IV, RV) RV = R CALL BYTPUT (B0, IV, 0) CALL BYTPUT (B1, IV, 1) CALL BYTPUT (B2, IV, 2) CALL BYTPUT (B3, IV, 3) CALL BYTGET (B3, IV, 0) CALL BYTGET (B2, IV, 1) CALL BYTGET (B1, IV, 2) CALL BYTGET (B0, IV, 3) R = RV RETURN END ! of swab32 *----------------------------------------------------------------------- * wpgm - write a Portable Gray Map * The maximum intensity is assumed to be 255. * *___Name______Type_____________In/Out___Description_____________________ * OFIL Character*(*) In File name * Y(LDY,M) Real In Image matrix * LDY Integer In Leading dimension of Y * N Integer In Columns * M Integer In Rows * IERR Integer Out Error code, 0 = no errors *----------------------------------------------------------------------- SUBROUTINE WPGM (OFIL, Y, LDY, N, M, IERR) ! write out IMPLICIT NONE INTEGER LDY, N, M, IERR REAL Y(LDY,M) CHARACTER*(*) OFIL INTEGER I, J, K, IOU, POS INTEGER*1 B, FID(3), MXV(4) CHARACTER*32 HEADER INTEGER IOUNIT, KTRIM, LTRIM ! external functions DATA FID / 80, 53, 10 / ! 'P', '5', Line feed DATA MXV / 50, 53, 53, 10 / ! '255', Line feed *----------------------------------------------------------------------- * Begin *----------------------------------------------------------------------- IOU = IOUNIT () OPEN (UNIT=IOU, FILE=OFIL, STATUS='NEW', ACCESS='DIRECT', RECL=1, $ FORM='UNFORMATTED', IOSTAT=IERR) IF (IERR .NE. 0) THEN PRINT *, 'trouble opening output file: ', OFIL RETURN END IF * Write 'P5' POS = 0 DO I = 1, 3 POS = POS + 1 WRITE (UNIT=IOU, REC=POS, IOSTAT=IERR) FID(I) IF (IERR .NE. 0) PRINT *, 'trouble writing ', FID(I) END DO * Time horizontal 10 FORMAT (I9) WRITE (UNIT=HEADER, FMT=10, IOSTAT=IERR) N IF (IERR .NE. 0) THEN PRINT *, 'trouble composing header!' RETURN END IF J = KTRIM(HEADER) K = LTRIM(HEADER) DO I = J, K POS = POS + 1 WRITE (UNIT=IOU, REC=POS, IOSTAT=IERR) HEADER(I:I) IF (IERR .NE. 0) PRINT *, 'trouble writing ', I, 'th', $ ' position of header' END DO * A single whitespace character B = 32 ! space POS = POS + 1 WRITE (UNIT=IOU, REC=POS, IOSTAT=IERR) B IF (IERR .NE. 0) PRINT *, 'trouble writing space!' * Frequency vertical WRITE (UNIT=HEADER, FMT=10, IOSTAT=IERR) M IF (IERR .NE. 0) THEN PRINT *, 'trouble composing header!' RETURN END IF J = KTRIM(HEADER) K = LTRIM(HEADER) DO I = J, K POS = POS + 1 WRITE (UNIT=IOU, REC=POS, IOSTAT=IERR) HEADER(I:I) IF (IERR .NE. 0) PRINT *, 'trouble writing ', I, 'th', $ ' position of header' END DO * A single whitespace character B = 10 ! line feed POS = POS + 1 WRITE (UNIT=IOU, REC=POS, IOSTAT=IERR) B IF (IERR .NE. 0) PRINT *, 'trouble writing line feed!' * Write '255' DO I = 1, 4 POS = POS + 1 WRITE (UNIT=IOU, REC=POS, IOSTAT=IERR) MXV(I) IF (IERR .NE. 0) PRINT *, 'trouble writing ', MXV(I) END DO * The data DO K = M, 1, -1 DO J = 1, N I = NINT(Y(J,K)) I = MIN(MAX(0, I), 255) CALL BYTPUT (B, I, 0) POS = POS + 1 WRITE (UNIT=IOU, REC=POS, IOSTAT=IERR) B IF (IERR .NE. 0) PRINT*, 'trouble writing bitmap at ', J, K END DO END DO * Done CLOSE (UNIT=IOU, IOSTAT=IERR) IF (IERR .NE. 0) PRINT *, 'trouble closing output file!' RETURN END ! of wpgm *----------------------- End of file spectr.f --------------------------