convolve_isa.h

convolve_isa.c

/*
 *  convolve_isa.c
 *
 *
 *  Created by Kenneth Buck on 12-11-27.
 *  Copyright 2012 University of Calgary. All rights reserved.
 *
 */

/* Include Files */
#include "convolve_isa.h"
#include "stdio.h"
#include "stdlib.h"
#include "math.h"

/* Function Prototypes */
void convolve(short x[], int N, float h[], int M, float y[], int P);
void print_vector(char *title, float x[], int N);

/************************************************************************
*
*	Function: main
*
*	Description: Main routine for the convolution of two files.
*
*************************************************************************/

int main(int argc, char *argv[]) {

	/* Ensure correct number of arguments */
	if(argc != 3) {
		printf("usage <inputFile> <IRfile> <outputFile>\n");
		return -1;
	}

	/***
	 * Input File
	 */
	FILE *inputFile = fopen( argv[1], "r" );
	struct waveHeader inputHeader;

	if(inputFile != 0) {
		fread(&inputHeader.chunkid,1,4,inputFile);
		fread(&inputHeader.chunksize,1,4,inputFile);
		fread(&inputHeader.format,1,4,inputFile);
		fread(&inputHeader.subchunk1id,1,4,inputFile);
		fread(&inputHeader.subchunksize,1,4,inputFile);
		fread(&inputHeader.audioformat,1,2,inputFile);
		fread(&inputHeader.numchannels,1,2,inputFile);
		fread(&inputHeader.samplerate,4,1,inputFile);
		fread(&inputHeader.byterate,4,1,inputFile);
		fread(&inputHeader.blockalign,2,1,inputFile);
		fread(&inputHeader.bitspersample,2,1,inputFile);
		fread(&inputHeader.subchunk2ID,4,1,inputFile);
		fread(&inputHeader.subchunk2size,4,1,inputFile);
	}

	// read actual data into an array.
	char inputData[inputHeader.subchunk2size];
	fread(&inputData,inputHeader.subchunk2size,1,inputFile);
	fclose( inputFile );

	// write data to a file

	int i;
	short val;

	// prepare input for 16 bit sample size.
	int inputSignalSize = inputHeader.subchunk2size / 2;
	short inputSignal[inputSignalSize];
	for(i=0;i<inputSignalSize;i+=2) {
		val = (short)((unsigned char) inputData[i]);
		val += (short)((unsigned char) inputData[i+1]) * 256;
		inputSignal[i/2] = val;
	}

	// normalize input
	float normalInputSignal[inputSignalSize];
	for(i = 0; i < inputSignalSize; i++ ) {
		normalInputSignal[i] = inputSignal[i] / pow(2,16);
	}

	/***
	 * Impulse File
	 */
	FILE *impulseResponse = fopen(argv[2], "r");
	struct waveHeader impulseHeader;

	if(impulseResponse != 0) {
		fread(&impulseHeader.chunkid,1,4,impulseResponse);
		fread(&impulseHeader.chunksize,1,4,impulseResponse);
		fread(&impulseHeader.format,1,4,impulseResponse);
		fread(&impulseHeader.subchunk1id,1,4,impulseResponse);
		fread(&impulseHeader.subchunksize,1,4,impulseResponse);
		fread(&impulseHeader.audioformat,1,2,impulseResponse);
		fread(&impulseHeader.numchannels,1,2,impulseResponse);
		fread(&impulseHeader.samplerate,4,1,impulseResponse);
		fread(&impulseHeader.byterate,4,1,impulseResponse);
		fread(&impulseHeader.blockalign,2,1,impulseResponse);
		fread(&impulseHeader.bitspersample,2,1,impulseResponse);
		fread(&impulseHeader.subchunk2ID,4,1,impulseResponse);
		fread(&impulseHeader.subchunk2size,4,1,impulseResponse);
	}

	// read impulse data into array.
	float impulseData[impulseHeader.subchunk2size];
	fread(&impulseData,impulseHeader.subchunk2size,1,impulseResponse);
	fclose( impulseResponse );

	// prepare input for 16 bit sample size.
	int impulseSignalSize = impulseHeader.subchunk2size / 2;
	short impulseSignal[impulseSignalSize];
	for(i=0;i<impulseSignalSize;i+=2) {
		val = (short)((unsigned char) impulseData[i]);
		val += (short)((unsigned char) impulseData[i+1]) * 256;
		impulseSignal[i/2] = val;
	}

	// normalize impulse
	float normalImpulseSignal[impulseSignalSize];
	for(i = 0; i < impulseSignalSize; i++ ) {
		normalImpulseSignal[i] = impulseSignal[i] / pow(2,16);
	}

	// create output array.(divide by two due to 16bits per sample)
	int outputSignalSize = inputSignalSize + impulseSignalSize - 1;
	float outputSignal[outputSignalSize];

	// do the convolution
	convolve(inputSignal,inputSignalSize,normalImpulseSignal,impulseSignalSize,outputSignal,outputSignalSize);

	// print result
	print_vector("output",outputSignal,outputSignalSize);

	return 0;
}

/************************************************************************
 *
 *	Function:	 convolve
 *
 *	Description: Convolves two signals, producing an output signal
 *				 The convolution is done in the time domain using
 *				 "Input Side Algorithm" (see Smitth, p. 112-115)
 *
 *	Parameters:	 x[] is the signal to be convolved
 *				 N is the number of samples in the vector x[]
 *				 h[] is the impulse response, which is convolved with x[]
 *				 M is the number of samples in the vector h[]
 *				 y[] is the output signal, the result of the convolution
 *				 P is the number of samples in the vector y[]. P must
 *					equal N + M - 1
 *
 *************************************************************************/

void convolve(short x[], int N, float h[], int M, float y[], int P)
{
	int n, m;

	/* Make sure the output buffer is the right size: P = N + M - 1 */
	if (P != (N + M - 1)) {
		printf("Output signal vector is the wrong size\n");
		printf("It is %-d, but should be %-d\n", P, (N + M - 1));
		printf("Aborting convolution\n");
		return;
	}

	/* Clear the output buffer y[] to all zero values */
	for (n = 0; n < P; n++)
		y[n] = 0.0;

	/* Do the convolution */
	/* Output loop: process each input value x[n] in turn */
	for (n = 0; n < N; n++){
		/* Inner loop: process x[n] with each sample of h[] */
		for (m = 0; m < M; m++) {
			y[n+m] += x[n] * h[m];
		}
	}
}

/************************************************************************
 *
 *	Function:	 print_vector
 *
 *	Description: Prints the vector out to the screen
 *
 *	Parameters:	 title is a string naming the vector
 *				 x[] is the vector to be printed out
 *				 N is the number of samples in the vector x[]
 *
 *************************************************************************/

void print_vector(char * title, float x[], int N)
{
	int i;

	printf("\n%s\n", title);
	printf("Vector size: %-d\n", N);
	printf("Sample Number \tSample Value\n");
	for(i = 0; i< N; i++)
		printf("%-d\t\t%f\n", i, x[i]);
}

gistfile1.c

/*
 * convolve_isa.h
 *
 *  Created on: Nov 29, 2012
 *      Author: kennethbuck
 */

#ifndef CONVOLVE_ISA_H_
#define CONVOLVE_ISA_H_

struct waveHeader {
	char 	chunkid[4];
	int 	chunksize;
	char 	format[4];
	char 	subchunk1id[4];
	int		subchunksize;
	short 	audioformat;
	short	numchannels;
	int		samplerate;
	int		byterate;
	short	blockalign;
	short 	bitspersample;
	char	subchunk2ID[4];
	int 	subchunk2size;
};

#endif /* CONVOLVE_ISA_H_ */