convolution and arithmetic mean

 
Hello,

I have a discrete signal a(x[i]) = cos(2*pi*x[i]/L).
I integrate analytically this cosine function from x[i-1] to x[i+1], I get

exact[i] =  L * (sin(2.0*pi*x(i+1)/L) - sin(2*pi*x(i-1)/L)) / ( 2.0*pi).

Now, I want to calculate the convolution a * h where h is a rectangular
filter (non causal filter). Doing that, I should recover for each node [i]
the exact value of the integral (exact[i]). 

I have attached my source code here below (Matlab code). I am doing
something wrong because I can get the ratio exact solution /numerical
solution to be equal to one. However,  this ratio is equal for each node
[i] to a constant value.

Does anyone have a idea ?

thank you,

manu


function fs=conv()

clc;
close all;
clear all;

L = 10;

n = 10;
m =  3;
p = n + m - 1;

x  = zeros(1,p);
dx = L/(p-2);
x  = [0:dx:L+dx];

% signal cos(x)
sigx      = cos(2.0*pi*x/L);
sigx(p)   = 0.0; //zero padding
sigx(p-1) = 0.0; //zero padding

% filter rect(x) //non causal filter
h    = zeros(1,p);
h(1) = 1.0/3.0;
h(2) = 1.0/3.0;
h(p) = 1.0/3.0;

%compute convolution in fourier Space
fx   = fft(sigx);
fh   = fft(h);
iconv = fx .* fh;
conv  = ifft(iconv);

%delete space used for zero-padding
conv(p)=[];
conv(p-1)=[];


% exact solution
exact = zeros(1,n);
xm1 = x - dx;
xp1 = x + dx;

for i=2:n-1,
 exact(i) =  L * (sin(2.0*pi*x(i+1)/L) - sin(2*pi*x(i-1)/L)) / ( 2.0*pi);
end


% display exact and numerical solution from convolution and ratio
exact
conv
ratio = exact ./ conv

On Wed, 9 Feb 2011 14:38:59 -0800 (PST), Clay <clay@claysturner.com>
wrote:

    [Snipped by Lyons]
>
>I'm not quite sure what you are trying to do. So I'll have to make
>some assumptions. First you want to integrate a discrete signal? Thus
>you need to "connect the dots" to make this into something you can
>integrate. If you connect the dots (nodes) with sinc functions, then
>the integrations becomes pretty simple since the integral of sinc (-
>infty to infty) is pi. And your infinite duration signal can now be
>integrated be realizing it is expressable as a sum of weighted sincs.
>So your integral is now a sum.
>
> Of course you have to scale your sinc functions so the zeroes of the
>sincs occur at the nodes (sample points). So this will result in a
>common scale factor for your whole integral.
>
>Now the hard part is what do you do if you are doing a definite
>integral. You can't just use the nodes inside of the interval since
>the sinc interpolation includes effects from the samples outside of
>the interval. I.e., you will need the integrals of sinc functions over
>parts of the domain. The end effects can be troublesome.
>
>I hope this helps. I suggest 1st trying a function that decays nicely
>to zero at both ends so then you can easily approximate the integral
>over the infinite domain by just using a large finite domain that
>covers most of the "energy" of the function to be integrated. Thus you
>minimize the end effects.
>
>
>Clay

Hi Clay,
   manuge2 is not integrating in any classical sense as far 
as I can tell.  He's *NOT* using Simpson's rule, Trapezoidal, 
nor Tick's rule for integration.  What he's doing 
is convolving his input sequence with the following 
sequence, "h":

    h = 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1

He does that by multiplying the DFT of his input 
sequence (padded with two zero-valued samples) 
times the DFT of the above "h" sequence, then performing an 
inverse DFT on that product.  So manuge2 is performing some 
sort of lowpass filtering (I think).  I'd have to study 
the behavior of that oddball "h" sequence before I could 
say any more about what manuge2 is actually doing.
In any case, I don't believe he's integrating his input.

See Ya',
[-Rick-]

>On Wed, 9 Feb 2011 14:38:59 -0800 (PST), Clay <clay@claysturner.com>
>wrote:
>
>    [Snipped by Lyons]
>>
>>I'm not quite sure what you are trying to do. So I'll have to make
>>some assumptions. First you want to integrate a discrete signal? Thus
>>you need to "connect the dots" to make this into something you can
>>integrate. If you connect the dots (nodes) with sinc functions, then
>>the integrations becomes pretty simple since the integral of sinc (-
>>infty to infty) is pi. And your infinite duration signal can now be
>>integrated be realizing it is expressable as a sum of weighted sincs.
>>So your integral is now a sum.
>>
>> Of course you have to scale your sinc functions so the zeroes of the
>>sincs occur at the nodes (sample points). So this will result in a
>>common scale factor for your whole integral.
>>
>>Now the hard part is what do you do if you are doing a definite
>>integral. You can't just use the nodes inside of the interval since
>>the sinc interpolation includes effects from the samples outside of
>>the interval. I.e., you will need the integrals of sinc functions over
>>parts of the domain. The end effects can be troublesome.
>>
>>I hope this helps. I suggest 1st trying a function that decays nicely
>>to zero at both ends so then you can easily approximate the integral
>>over the infinite domain by just using a large finite domain that
>>covers most of the "energy" of the function to be integrated. Thus you
>>minimize the end effects.
>>
>>
>>Clay
>
>Hi Clay,
>   manuge2 is not integrating in any classical sense as far 
>as I can tell.  He's *NOT* using Simpson's rule, Trapezoidal, 
>nor Tick's rule for integration.  What he's doing 
>is convolving his input sequence with the following 
>sequence, "h":
>
>    h = 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1
>
>He does that by multiplying the DFT of his input 
>sequence (padded with two zero-valued samples) 
>times the DFT of the above "h" sequence, then performing an 
>inverse DFT on that product.  So manuge2 is performing some 
>sort of lowpass filtering (I think).  I'd have to study 
>the behavior of that oddball "h" sequence before I could 
>say any more about what manuge2 is actually doing.
>In any case, I don't believe he's integrating his input.
>
>See Ya',
>[-Rick-]
>

Hello,

Thank you for your comments and advice. Actually, I have replaced the
coefficients in the Kernel function so that they now correspond to the
Simpson's coefficients : i.e,  h = [4/6 1/6 0 0 0 0 0 0 0 0 0 0 1/6]

Using Simpson's method for integrating is equivalent to what I have done.

Thanks again,

Manu