rlm@23: #include rlm@23: #include rlm@23: #include rlm@23: #include rlm@23: #include rlm@23: #include "SndfileWavUtil.h" rlm@23: rlm@23: const int Taps = 9; rlm@23: const int Points = 8; rlm@23: const int log_Points = 3; rlm@23: static int newest; rlm@23: static short int x[Taps] = {0,0,0,0,0,0,0,0,0}; rlm@23: static short int fir(short int); rlm@23: rlm@23: // foreward declaration. This is defined in DFT.cpp rlm@23: int DFT(int dir,int m,double *x1,double *y1); rlm@23: rlm@23: rlm@23: static int x_to_y(int x, int y) rlm@23: { rlm@23: assert(y>0); rlm@23: int rv = x; rlm@23: while(y-->1) rlm@23: x *= x; rlm@23: return rv; rlm@23: } rlm@23: rlm@23: // function implemented in LinearWindowingFunction.bsv rlm@23: static double windowing_fn(int index) rlm@23: { rlm@23: double i = (double)index; rlm@23: double divisor = Points/2; //(double) (x_to_y(2,log_Points-1)); rlm@23: double rv = (index < Points/2) ? (i / divisor) : ((((double)Points)-i) / divisor); rlm@23: return rv; rlm@23: } rlm@23: rlm@23: // lifted from the bluespec code rlm@23: static double h[Taps] = {-0.0124, rlm@23: 0.0, rlm@23: -0.0133, rlm@23: 0.0, rlm@23: 0.8181, rlm@23: 0.0, rlm@23: -0.0133, rlm@23: 0.0, rlm@23: -0.0124}; rlm@23: rlm@23: rlm@23: short int rlm@23: fir(short int sample) rlm@23: { rlm@23: x[newest] = sample; rlm@23: double y = 0; rlm@23: for (int k = 0; k < Taps; k++) { rlm@23: short int bits = x[(newest+k)%Taps]; rlm@23: double x = (double)bits; rlm@23: y += h[k] * x; rlm@23: } rlm@23: newest = (newest == 0) ? (Taps-1) : (newest-1); rlm@23: short int rv = ((short int)y)&0x0000FFFF; rlm@23: return rv; rlm@23: } rlm@23: rlm@23: rlm@23: int rlm@23: main (int argc, char * argv []) rlm@23: { rlm@23: const char* inputWavFileName; rlm@23: const char* outputWavFileName; rlm@23: rlm@23: FILE *inputPcmFile; rlm@23: FILE *firPcmFile; rlm@23: FILE *ifftPcmFile; rlm@23: FILE *outputPcmFile; rlm@23: rlm@23: short int sample; rlm@23: short int samples[Points]; rlm@23: const unsigned int short_sz = sizeof(short int); rlm@23: rlm@23: inputWavFileName = argv[1]; rlm@23: outputWavFileName = argv[2]; rlm@23: rlm@23: // Convert input wav to pcm rlm@23: generate_pcm(inputWavFileName, "input.pcm"); rlm@23: rlm@23: inputPcmFile = fopen("input.pcm", "r"); rlm@23: firPcmFile = fopen("fir.pcm", "w"); rlm@23: rlm@23: assert(inputPcmFile); rlm@23: assert(firPcmFile); rlm@23: rlm@23: newest = 0; rlm@23: memset(x, 0, sizeof(x)); rlm@23: rlm@23: while(fread(&sample, short_sz, 1, inputPcmFile)) { rlm@23: sample = fir(sample); rlm@23: assert(fwrite(&sample,short_sz, 1, firPcmFile)); rlm@23: } rlm@23: rlm@23: fclose(inputPcmFile); rlm@23: fclose(firPcmFile); rlm@23: rlm@23: firPcmFile = fopen("fir.pcm", "r"); rlm@23: ifftPcmFile = fopen("ifft.pcm", "w"); rlm@23: rlm@23: assert(firPcmFile); rlm@23: assert(ifftPcmFile); rlm@23: rlm@23: int read = 0; rlm@23: rlm@23: // read in half a frame rlm@23: for(int i = 0; i < Points/2; i++) rlm@23: read += fread(&samples[(Points/2)+i], short_sz, 1, firPcmFile); rlm@23: rlm@23: if(read > 0) rlm@23: assert(read==Points/2); rlm@23: rlm@23: // we will do an 'Points' point fft, and then undo it rlm@23: while(true) { rlm@23: rlm@23: read = 0; rlm@23: rlm@23: // shift last half frame to head rlm@23: for(int i = 0; i < Points/2; i++) rlm@23: samples[i] = samples[(Points/2)+i]; rlm@23: rlm@23: // read in another half frame rlm@23: for(int i = 0; i < Points/2; i++) rlm@23: read += fread(&samples[(Points/2)+i], short_sz, 1, firPcmFile); rlm@23: rlm@23: if(read == 0) rlm@23: break; rlm@23: rlm@23: // pad out the rlm@23: if(read < Points/2){ rlm@23: for(int i = read; i < Points/2; i++) rlm@23: samples[Points/2+i] = 0; rlm@23: } rlm@23: rlm@23: double dsamples[Points]; rlm@23: double dimag[Points]; rlm@23: rlm@23: // this shift is performed in the Bluespec rlm@23: for(int i = 0; i < Points; i++){ rlm@23: dsamples[i] = (double)(samples[i]>>log_Points); rlm@23: dimag[i] = 0.0; rlm@23: } rlm@23: rlm@23: rlm@23: DFT(1,Points,dsamples,dimag); rlm@23: DFT(-1,Points,dsamples,dimag); rlm@23: rlm@23: short int ifftResult[Points]; rlm@23: for(int i = 0; i < Points; i++) rlm@23: ifftResult[i] = (int)dsamples[i]; //)>>log_Points; rlm@23: rlm@23: // strip off the padding rlm@23: int write = (read < Points/2) ? read+(Points/2) : Points; rlm@23: for(int i = 0; i < write; i++) rlm@23: assert(fwrite(&ifftResult[i],short_sz,1,ifftPcmFile)); rlm@23: rlm@23: // break out if we're at the end rlm@23: if(read < Points/2){ rlm@23: break; rlm@23: } rlm@23: } rlm@23: rlm@23: fclose(firPcmFile); rlm@23: fclose(ifftPcmFile); rlm@23: rlm@23: // do the windowing rlm@23: rlm@23: ifftPcmFile = fopen("ifft.pcm", "r"); rlm@23: outputPcmFile = fopen("output.pcm", "w"); rlm@23: rlm@23: assert(ifftPcmFile); rlm@23: assert(outputPcmFile); rlm@23: rlm@23: short int samplesA[Points]; rlm@23: short int samplesB[Points]; rlm@23: rlm@23: bool valid_stream = true; rlm@23: // write out the first half frame rlm@23: for(int i = 0; i < Points/2; i++){ rlm@23: valid_stream &= fread(samplesA, short_sz, 1,ifftPcmFile); rlm@23: valid_stream &= fwrite(samplesA, short_sz,1,outputPcmFile); rlm@23: } rlm@23: rlm@23: while(valid_stream){ rlm@23: rlm@23: for(int i = 0; i < Points/2; i++) rlm@23: valid_stream &= fread(&samplesA[i], short_sz, 1,ifftPcmFile); rlm@23: rlm@23: for(int i = 0; i < Points; i++) rlm@23: valid_stream &= fread(&samplesB[i], short_sz, 1,ifftPcmFile); rlm@23: rlm@23: for(int i = 0; i < Points/2; i++) rlm@23: valid_stream &= fread(&samplesA[(Points/2)+i], short_sz, 1,ifftPcmFile); rlm@23: rlm@23: rlm@23: // this isn't quite right rlm@23: if(!valid_stream) rlm@23: break; rlm@23: rlm@23: for(int i = 0; i < Points; i++){ rlm@23: double window = windowing_fn(i); rlm@23: short int rv = samplesA[i]*(1-window) + window*samplesB[i]; rlm@23: assert(fwrite(&rv,short_sz,1,outputPcmFile)); rlm@23: } rlm@23: rlm@23: } rlm@23: rlm@23: fclose(ifftPcmFile); rlm@23: fclose(outputPcmFile); rlm@23: rlm@23: generate_wav("output.pcm", inputWavFileName, outputWavFileName); rlm@23: rlm@23: } rlm@23: