BobM wrote:
> It's difficult what to say what is wrong without knowing more
> information about how your function works. What is the total input size
> in samples? The filter size? The input block size?
I used two test input signals, one 2,103,154 samples long, the other,
972 samples long. as for the filter, it's 20,000 samples long in the
case of the shifted delta function. I'm not sure to know what you mean
by input block size
> I'm assuming your filter size is 20,000 taps. If my math is correct,
> the max your input block size can be is 12,769 samples (32768-20000+1).
huh? ok, we maybe got the finger on something wrong... what I do is, in
the case of the convolution by the shifted delta function, take the
20,000 samples long signal, and provided that the other input signal is
longer than that, zero-pad that signal to the power of two above twice
the length of the signal (I guess that is 65,536 here), then, I cut the
2,103,154 samples long signal into blocks of 32,768 samples and I
zero-pad each block to 65,536 (i think i made no mistake)
> Anyway, the result you're getting seems correct to me. As long as your
> input signal is shifted in time by 4000 samples, which it sounds like
> it is, then this is the correct behavior.
No, it's not getting shifted, but muted, I mean, some regions of 4000
samples are getting, each 32,768 samples, set to zero, as the other
28,768 samples are unchanged, even when using a windowed-sinc function.
This is indeed quite far from a normal behavior
> Also, it seems your input
> block isn't too big, because I'd expect to see some non-zero samples
> folding back into the first 4000 samples if it was too big. Sounds like
> the convolution itself is working.
>
> > I also tried with a windowed-sinc function, i used 88.2 Hz as a cutoff
> > frequency, 8.82 Hz as the bandwidth of the rolloff, and when convolving
> > with the input signal, not only the output signal has regions regularly
> > set to zero, but the output signal doesnt look/sound at all like it's
> > been low-passed filtered
>
> How many taps is this filter?
20,001. I calculate the length of the filter by 4/BW (BW=the bandwidth
of the roll off) and make sure the number is odd.
> If it's not many taps, then this may be
> such a shallow filter that it there's almost no effect.
> What kind of
> window is being used to generate the filter?
Blackman. I made sure it's working right. I even made sure that the
first and last bins were always non-zero.
> According to the first
> result you mentioned it sounds like your convolution is working, so
> this may be related to the filter you're using. Does this same filter
> work with a time-domain implementation of convolution?
Haven't made a time-domain convolution. i thought i'd go straight to
the goal... anyways, it looks perfectly right.
>
> > any ideas on what can be wrong?
>
> Again, please provide some more information. However, my hunch is that
> the convolution itself is working. But confirm that your input block
> isn't too big. Also check the way you're piecing the output blocks back
> together. Whether you're doing overlap-add or overlap-save, check that
> the overlapped portions are being added together before sending them
> off to output.
Maybe it's not very appropriate to paste such a long code, but I have
the feeling that it'll go much faster if you can get to see a
stripped-down-to-the-essential version of my function, so here it is
(the first functions are here because they are used in the fftconv
function)
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <math.h>
#include <fftw3.h>
double log_2(double x) //log base 2 performing function
{
double y;
y=log(x)/log(2);
return y;
}
double modfloat(double a, double b) //modulo returning function
{
double c, d, x;
c=(uint32_t) a/b;
d=a/b;
x=d-c;
return x;
}
double roundup(double in) //function rounding the input to the upper
integer (unless the input already is an integer)
{
double out;
if (modfloat(in,1) == 0)
out=in;
else
out=in - (modfloat(in,1)) + 1;
return out;
}
double *zeropad(double *s, uint32_t nin, uint32_t nout)
{
uint32_t i;
double *h;
h=malloc (sizeof(double) * nout);
for (i=0; i<nin; i++)
h[i]=s[i];
for (i=nin; i<nout; i++)
h[i]=0;
return h;
}
void dftfunction(double *in, uint32_t N, uint8_t method, uint8_t root)
{
uint32_t i;
fftw_plan p;
p = fftw_plan_r2r_1d(N, in, in, method, FFTW_ESTIMATE);
fftw_execute(p);
fftw_destroy_plan(p);
if (root==1)
for (i=0;i<N;i++)
in[i]=in[i]/sqrt(N);
}
double *fftconv(double *s1, uint32_t n1, double *s2, uint32_t n2,
uint32_t *n3)
//s1 and s2 the two input signals, n1 and n2 their respective length
{
uint32_t i, j, nmin, nmax, blocksize, nblocks;
double *s3, *smin, *smax, **block, a, b, c, d;
printf("fftconv...\n");
*n3=n1+n2-1; //output signal length
s3=malloc (*n3 * sizeof(double));
if (n1<n2)
{
smin=s1;
smax=s2;
nmin=n1;
nmax=n2;
}
else
{
smin=s2;
smax=s1;
nmin=n2;
nmax=n1;
}
blocksize=pow(2, (uint32_t) roundup(log_2((double) (nmin)*2))); //size
of blocks for the overlap-add of the long signal, size of the
zeropadded short signal
smin=zeropad(smin, nmin, blocksize);
dftfunction(smin, blocksize, 0, 0); //FFT of the shortest signal
nblocks=*n3/(blocksize/2); //number of blocks for the o-a of the long
signal
block=malloc (nblocks * sizeof(double *));
for (i=0; i<nblocks; i++) //loop filling the contents of the long
signal in an array of blocks
{
block[i]=malloc (blocksize * sizeof(double));
for (j=0; j<blocksize/2; j++)
block[i][j]=smax[i*(blocksize/2) + j];
for (j=blocksize/2; j<blocksize; j++)
block[i][j]=0;
dftfunction(block[i], blocksize, 0, 1); //FFT of each block
}
for (i=0; i<nblocks; i++)
{
a=block[i][0]; // DC
c=smin[0];
block[i][0]=a*c; // no need for imaginary parts
for (j=1; j<blocksize/2; j++)
{
a=block[i][j];
b=block[i][blocksize-j];
c=smin[j];
d=smin[blocksize-j];
block[i][j]=a*c - b*d;
block[i][blocksize-j]=b*c + a*d;
}
a=block[i][blocksize/2]; // Nyquist
c=smin[blocksize/2];
block[i][blocksize/2]=a*c;
dftfunction(block[i], blocksize, 1, 1); //IFFT of the result of the
multiplication in each block
}
for (i=0; i<*n3; i++) //initialisation of the output signal s3
s3[i]=0;
for (i=0; i<nblocks; i++) //filling of s3 from the content of the
blocks
for (j=0; j<blocksize/2; j++)
s3[i*(blocksize/2) + j]=s3[i*(blocksize/2) + j]+block[i][j];
return s3;
free(block);
}