Adaptive Filter | ADPCM | ADSP | ADSP-2181 | Aliasing | AMR | Anti-Aliasing | ARMA | Autocorrelation | AutoCovariance | Beamforming | Bessel | Blackfin | Butterworth | C6713 | CCS | Chebyshev | CIC Filter | Circular Convolution | Code Composer Studio | Comb Filter | Compression | Convolution | Cross Correlation | DCT | Decimation | Deconvolution | Demodulation | DM642 | DSP Boards | DSP/BIOS | DTMF | Echo Cancellation | Equalization | Equalizer | ETSI | EZLITE (Ez-kit Lite) | FFT | FFTW | FIR Filter | Fixed Point | FSK | G.711 | G.723 | G.729 | Gaussian Noise | Goertzel | GPIO | Hilbert Transform | IFFT | IIR Filter | Interpolation | Invariance | JTAG | Kalman | Laplace Transform | Levinson | LPC | McBSP | MIPS | Modulation | MPEG | Multirate | Notch Filter | Nyquist | OFDM | Oversampling | Pink Noise | Pitch | PLL | Polyphase | QAM | QDMA | Quantization | Quantizer | Radar | Random Noise | Reed Solomon | Remez | Resampling | RTDX | Sampling | Sharc | TI C6711 | Undersampling | Viterbi | Wavelets | White Noise | Wiener Filter | Windowing | XDS510PP | Z Transform

So, here is my problem.

(I have put some matrices here, because line breaks in usenet, would
make it hard to read: www.fledelius.com/data.html)

I take a gaussian function in the spatial domain with i=0.
I take an input with one entry=1 and the rest 0 and i=0.
I convert both to the fourier domain and multiply them.
I take the inverse FFT of the result.

In theory this should be the equvalent of convolving with a gaussian
kernel, right?

If you look at the last matrix at the bottom of the page, it looks like
the result in spatial domain have swapped quadrants. 

I could of course unswap them, but I need to this in 3D, so not sure it
is the right way to go?

Any thoughts, thank you,
Walther Fledelius

(I use FFTW3, XCode, OS X 10.3.2)

______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.

Walther Fledelius wrote:
> So, here is my problem.
> 
> (I have put some matrices here, because line breaks in usenet, would
> make it hard to read: www.fledelius.com/data.html)
> 
> I take a gaussian function in the spatial domain with i=0.
> I take an input with one entry=1 and the rest 0 and i=0.
> I convert both to the fourier domain and multiply them.
> I take the inverse FFT of the result.
> 
> In theory this should be the equvalent of convolving with a gaussian
> kernel, right?

Walther,

You're not listening. Multiplying using the FFT (which is just an
efficient DFT) is performing *circular* convolution. Because of
that, you have to ensure certain properties of the input sequences
are met (like padding with a certain number of zeros) to ensure
that this gives you plain old convolution. Otherwise you get "time
aliasing." Read P&M like I suggested.
-- 
%  Randy Yates                  % "...the answer lies within your soul
%% Fuquay-Varina, NC            %       'cause no one knows which side
%%% 919-577-9882                %                   the coin will fall."
%%%% <y...@ieee.org>           %  'Big Wheels', *Out of the Blue*, ELO
http://home.earthlink.net/~yatescr

______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.

In article <CnYFb.6011$I...@newsread3.news.atl.earthlink.net>,
Randy Yates <y...@ieee.org> wrote:

> Walther Fledelius wrote:
> > So, here is my problem.
> > 
> > (I have put some matrices here, because line breaks in usenet, would
> > make it hard to read: www.fledelius.com/data.html)
> > 
> > I take a gaussian function in the spatial domain with i=0.
> > I take an input with one entry=1 and the rest 0 and i=0.
> > I convert both to the fourier domain and multiply them.
> > I take the inverse FFT of the result.
> > 
> > In theory this should be the equvalent of convolving with a gaussian
> > kernel, right?
> 
> Walther,
> 
> You're not listening. Multiplying using the FFT (which is just an
> efficient DFT) is performing *circular* convolution. Because of
> that, you have to ensure certain properties of the input sequences
> are met (like padding with a certain number of zeros) to ensure
> that this gives you plain old convolution. Otherwise you get "time
> aliasing." Read P&M like I suggested.

Oh, but I am, and I have been for a few days, reading a lot as well!
My bookstore is closed for a few days, so bare with me.

I will do google, now that I know the circularity might be the problem,
maybe I can find the answer there.

Thank you,
Walther

______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.

Walther Fledelius wrote:
...
> My bookstore is closed for a few days, so bare with me.

Try this (Chapter 18 - FFT Convolution):

http://www.dspguide.com/pdfbook.htm

Regards,
Andor

______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.

"Walther Fledelius" <-...@-.dk> wrote in message
news:231220031625598959%-...@-.dk...
> In article <CnYFb.6011$I...@newsread3.news.atl.earthlink.net>,
> Randy Yates <y...@ieee.org> wrote:
>
> > Walther,
> >
> > You're not listening. Multiplying using the FFT (which is just an
> > efficient DFT) is performing *circular* convolution. Because of
> > that, you have to ensure certain properties of the input sequences
> > are met (like padding with a certain number of zeros) to ensure
> > that this gives you plain old convolution. Otherwise you get "time
> > aliasing." Read P&M like I suggested.
>
> Oh, but I am, and I have been for a few days, reading a lot as well!
> My bookstore is closed for a few days, so bare with me.
>
> I will do google, now that I know the circularity might be the problem,
> maybe I can find the answer there.

Walther,

A couple of suggestions:

1) In reviewing your data tables: Complex multiplication yields a complex
result even if one of the multiplicands has zero imaginary part.  So, you're
going to have to deal with that.
p=(a+jb)x(c+jd)= (ac-bd)+j*(bc+ad) so if d=0, p still has an imaginary term
jbc; p=ac+jbc
So, the final result you get in the frequency domain will not have zero
imaginary part.

2) Circular convolution requires that you zero-pad the data before taking
the Fourier Transform.  You've not done that.  The dimensions need to be at
least M+N-1 for no overlap.  So an MxM matrix would become (M+N-1,M+N-1) and
the same for an NxN matrix.  For M=2 and N=3:
M= a b   N=e f g
   c d     h i j
           k l m

However, a circular convolution would yield a 3x3, not a 4x4

I'm going to expand N to a repeating 4x4 just so I don't lose my mind in
doing the convolution:

Nx= e f g e
    h i j h
    k l m k
    e f g e

Now, the circular convolution of M and N can be computed using only the full
overlaps of M on Nx

ae+bf+ch+di  af+bg+ci+dj  ag+be+cj+dh
ah+bi+ck+dl  ai+bj+cl+dm  aj+bh+cm+dk
ak+bl+ce+df  al+bm+cf+dg  am+bk+cg+de

Notice that there are term here that wrap around.

Now, let's look at a circular convolution that avoids aliasing or wrap
around.  First, we zero-pad the matrices to 4x4 - which is the dimension
needed for no overlap:

M'=0 0 0 0  a 4x4
   0 a b 0
   0 c d 0
   0 0 0 0
N'=e f g 0  a 4x4
   h i j 0
   k l m 0
   0 0 0 0

Again, to maintain my sanity, I will repeat N' so I can do this in my head:

N'x=e f g 0 e f g
    h i j 0 h i j
    k l m 0 k l m
    0 0 0 0 0 0 0
    e f g 0 e f g
    h i j 0 h i j
    0 0 0 0 0 0 0
Convolving by M' with full overlap and unless I made a mistake, yields:

ai+bj+cl+dm  aj+ 0+cm+ 0   0+bh+ 0+dk  ah+bi+ck+dl
al+bm+ 0+ 0  am+ 0+ 0+ 0   0+bk+ 0+ 0  ak+bl+ 0+ 0
 0+ 0+cf+dg   0+ 0+cg+ 0   0+ 0+ 0+de   0+ 0+ce+df
af+bg+ci+dj  ag+ 0+cj+ 0   0+be+ 0+dh  ae+bf+ch+di
(I've not shown the zero terms from M', only the zero terms from N')

So, this result has 7 more unique, nonzero results in it.  Also, there are
no "overlapped" terms such as
am+bk+cg+de.

Doing this by multiplication in the frequency domain is no different
conceptually.  So, you need to transform arrays like M' and N' to begin
with.

Note that I have not done anything about absolute registration other than to
put the zeros where I put them.  The results will vary according to where
you put the spatial 0,0 data.

Also note, if you convolve two real arrays, the result is real.  So, if you
multiply their (complex) transforms then the inverse transform of their
frequency domain product must still be real (less small imaginary values due
to numerical inaccuracy).

I hope this helps.

Fred

______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.

> I hope this helps.
> 
> Fred

Is does, very nice explanation! 

Thank you all for the help, I will try and change my code as you
suggested.

Merry Christmas,
Walther Fledelius


(once again, usenet is the place to meet helpful people who know their
field)

______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.