Hi,
I'm trying to do a simple effect in the frequency domain, nothing extraordinary.
I need to apply various gains to specific frequencies.
For this purpose, I want to process a 16bit audio input stream, block by block,
apply a hanning window, do a fft, now apply my effect, and then an inverse fft,
followed by overlap and add buffering.
For now, I'm simply trying to do all of this without applying any effect, just
to see if my time to frequency and back conversion is correct.
I managed to do this in Octave. Now, I'm trying to do it in C with fftw and the
sound I obtain is like "shivering". It's recognizable but severely altered. I've
spent hours on this, it just won't work.
Attached :
- freqfilter.m: my working octave version
- transform.h and transform.c: time -> frequency -> time routines
- transformtest.c: test program
Could anyone tell me what's wrong with this C code?
Thanks in advance
--
Olivier
blocksize=2048;
args=argv();
if length(args) != 3
printf('Usage: infile nframes outfile\n');
printf('Remark: use nframes=-1 to read the full infile\n');
exit(1)
end
[signal, samplerate, bits] = wavread(args{1});
signal = signal(:,1);
signallen = str2num(args{2});
if signallen < 0 || signallen > length(signal)
signallen = length(signal);
end
printf('read signal of length %d\n', length(signal));
printf('processing %d frames\n', signallen);
maxlevel=max(abs(signal))
output=[];
window=hanning(blocksize);
for start = 1:blocksize/2:signallen
stop = start + blocksize - 1;
block = zeros(blocksize,1);
if stop > signallen
stop = signallen;
end
block(1:stop - start + 1) = signal(start:stop);
block = block .* window;
freq_block = fft(block, blocksize)(1:blocksize / 2);
% effect
% for ex: freq_block(64:blocksize/2) = 0;
block = real(ifft(freq_block, blocksize));
if length(output) < start + blocksize - 1
output(length(output) + 1: start + blocksize - 1) = 0;
end
output(start:start + blocksize - 1) += transpose(block);
end
output = transpose(output);
gain = maxlevel / max(abs(output));
printf('applying gain: %f\n', gain);
output = output * gain;
printf('writing output');
wavwrite(output, samplerate, bits, args{3});
#ifndef TRANSFORM_H
#define TRANSFORM_H
#include <stdint.h>
#include <fftw3.h>
typedef struct Transform Transform;
typedef void (* TransformCallback)(fftwf_complex *freqdata, void *data);
Transform * transform_new(int blocksize, TransformCallback callback, void *data);
void transform_destroy(Transform *self);
/* Convert blocksize input frames and calls callback with the freqdata */
void transform_push(Transform *self, int16_t *timedata);
/* Return blocksize/2 output frames */
int16_t * transform_pop(Transform *self);
#endif // TRANSFORM_H
#include <stdlib.h>
#include <math.h>
#include "transform.h"
struct Transform {
int blocksize;
TransformCallback callback;
void * callback_data;
fftwf_plan fft;
fftwf_plan ifft;
float * overlap;
int16_t * outframes;
int windex;
int rindex;
int rlen;
float * hanning;
float * real;
fftwf_complex * cplx;
int begun;
float scale;
};
static void generate_hanning(float *window, int size) {
int i;
float factor = 2.f * M_PI / (size - 1.f);
for (i = 0; i < size; i++) {
window[i] = (1.f - cos(i * factor)) / 2.f;
}
}
Transform *transform_new(int blocksize, TransformCallback callback, void *data) {
Transform *self = calloc(1, sizeof(Transform));
self->blocksize = blocksize;
self->callback = callback;
self->callback_data = data;
self->cplx = fftwf_malloc((blocksize / 2 + 1) * sizeof(*self->cplx));
self->real = fftwf_malloc(blocksize * sizeof(*self->real));
self->fft = fftwf_plan_dft_r2c_1d(blocksize, self->real, self->cplx, FFTW_ESTIMATE);
self->ifft = fftwf_plan_dft_c2r_1d(blocksize, self->cplx, self->real, FFTW_ESTIMATE);
self->overlap = calloc(blocksize + blocksize / 2, sizeof(*self->overlap));
self->outframes = calloc(blocksize / 2, sizeof(*self->outframes));
self->hanning = calloc(blocksize, sizeof(*self->hanning));
generate_hanning(self->hanning, blocksize);
self->scale = INT16_MAX / (float) blocksize;
return self;
}
void transform_destroy(Transform *self) {
fftwf_free(self->real);
fftwf_free(self->cplx);
free(self->hanning);
free(self->overlap);
fftwf_destroy_plan(self->ifft);
fftwf_destroy_plan(self->fft);
free(self);
}
void transform_push(Transform *self, int16_t *frames) {
int i, ii;
for (i = 0; i < self->blocksize; i++) {
self->real[i] = frames[i] / (float) INT16_MAX * self->hanning[i];
}
fftwf_execute(self->fft);
self->callback(self->cplx, self->callback_data);
fftwf_execute(self->ifft);
ii = self->blocksize / 2;
float *overlap = self->overlap + self->windex * self->blocksize / 2;
for (i = 0; i < ii; i++) {
overlap[i] += self->real[i] * self->scale;
}
if (self->windex == 2) {
overlap = self->overlap;
self->windex = 0;
} else {
overlap += ii;
self->windex++;
}
float *real2 = self->real + ii;
for (i = 0; i < ii; i++) {
overlap[i] = real2[i] * self->scale;
}
if (!self->begun) {
self->begun = 1;
} else {
self->rindex = (self->windex + 1) % 3;
self->rlen = 2;
}
}
int16_t * transform_pop(Transform *self) {
int16_t *outframes = NULL;
if (self->rlen) {
float *overlap = self->overlap + self->rindex++ * self->blocksize / 2;
self->rindex %= 3;
self->rlen--;
outframes = self->outframes;
int i;
for (i = 0; i < self->blocksize / 2; i++)
outframes[i] = overlap[i];
}
return outframes;
}
/*
* Transform object test
*
* Compile with:
* gcc -Wall -o transformtest transformtest.c transform.c -lsndfile -lfftw3f
*/
#include <stdio.h>
#include <sndfile.h>
#include <stdlib.h>
#include <string.h>
#include "transform.h"
#define BLOCKSIZE 2048
static void transformed(fftwf_complex *freqdata, void *data) {
}
int main(int argc, char *argv[]) {
if (argc != 3) {
printf("Usage: infile outfile\n");
return 1;
}
SF_INFO info;
info.format = 0;
SNDFILE *infile = sf_open(argv[1], SFM_READ, &info);
SNDFILE *outfile = sf_open(argv[2], SFM_WRITE, &info);
int16_t *buffer = calloc(BLOCKSIZE, sizeof(*buffer));
Transform *transform = transform_new(BLOCKSIZE, transformed, NULL);
while (1) {
int read = sf_read_short(infile, buffer, BLOCKSIZE);
if (read < BLOCKSIZE) {
memset(buffer + read, 0, (BLOCKSIZE - read) * sizeof(*buffer));
}
transform_push(transform, buffer);
int16_t *processed;
while ((processed = transform_pop(transform)) != NULL) {
sf_write_short(outfile, processed, BLOCKSIZE / 2);
}
if (read < BLOCKSIZE) {
// losing the last half block
break;
}
}
transform_destroy(transform);
sf_close(outfile);
sf_close(infile);
return 0;
}
_______________________________________________
Linux-audio-dev mailing list
[email protected]
http://lists.linuxaudio.org/mailman/listinfo/linux-audio-dev
