# Tag Info

## Hot answers tagged convolution

16

EDIT: graphical alignment issue corrected. I think this has promise. The centering is almost perfect, and while the code may not be well tuned or optimized it is quite fast: 0.008736 seconds on my machine. It works by attempting to find the center of each white "blob" and then averaging those positions. img = Import["http://i.stack.imgur.com/i050B.png"]; ...

16

Your image: img = Import["http://i.stack.imgur.com/i050B.png"]; Here is the circle's radius and center: crcl = ComponentMeasurements[RegionBinarize[Dilation[img, 3], {{320, 240}}, 0.3], {"Centroid", "EquivalentDiskRadius"}]; // AbsoluteTiming crcl {0.126007, Null} {1 -> {{284.448, 241.873}, 109.256}} And here how precise it is: ...

13

All the image samples you posted are basically lines pointing towards a common center. We can turn that into a simple mathematical model: Let's say every gradient in the image is normal to a line pointing towards the center. Find the center with the least squared error. So the error term would be: squaredError = ({cx - x, cy - y}.{gx, gy})^2; where cx/cy ...

12

Use the following: ListConvolve[a, b, {1, -1}, 0] (* ==> {1, 3, 6, 10, 10, 10, 9, 7, 4} *) This says: align the first element of b with the 1st element of a, align the last element of b with the -1st (i.e. last) element of a, pad with 0s if necessary. Have you tried searching the docs for "convolution"?

11

I hope I see the essence here. You are interested in the convolution of an interpolated function with a Gauss function Your underlying data has regular spacings in x-direction and the convolution with a Gaussian is extremely fast implemented in GaussianFilter for discrete data. Why are you making it so complicated when the only thing you have to do is ...

10

I think this should work: ClearAll[r]; r[0, t_] := Exp[-k*t]*Cos[t]; r[n_, t_] := Integrate[r[0, t - td]*r[n - 1, td], {td, 0, t}] eg r[2,t] (* (\[ExponentialE]^(-k t) (2 \[ExponentialE]^(k t) k^2 - k (2 k + t + k^2 t) Cos[t] + (k - k^3 + t + k^2 t) Sin[t]))/(2 (1 + k^2)^2) *)

10

First I want to say, as you mentioned in your comment that your ultimate goal is to to do it for nMax over 100, I suggest you first symbolicly calculate the correlation of the following function, treating $r_n$ ($n=-s,-s+1,\dots,s$, and $s$ is nSteps for short) as variables as $x$: $$\xi(x,r_{-s},r_{-s+1},\dots,r_{s})=\sum _{n=-s}^{s} r_n\, ... 9 Reading a bit about convolutions and correlations will help. The correlation of functions f(x) and g(x) is defined as$$c(t) = \int_{-\infty}^{\infty} f(t+x) g(x) \; dx ImageCorrelate and ListCorrelate do the same thing for lists of numbers. For example, ListCorrelate[{a, b, c}, {1, 2, 3, 4, 5, 6, 7, 8, 9, 10}] returns a list whose first element ...

8

Here are two things you can do to speed up this code. 1. Do the convolution with symbolic y Because you have defined corr using SetDelayed, the table of Convolve expressions will be re-evaluated every time you evaluate corr[number]. The normalisation term with y=0 is causing a particular slow down, though I'm not sure why exactly. If you instead use Set ...

8

As already mentioned, this is a convolution. Luckily, there's a more natural function to use for this problem than Integrate[], and that function is called, appropriately enough, Convolve[]. Now, since Convolve[] assumes an infinite integration region, we need a UnitStep[] multiplier in both the functions being convolved to limit the integration region to a ...

8

In line with the OPs request for a comparison of several different methods here's a comparison of six different ways to filter (or convolve) a data set x with a kernel h: convolution, correlation, the frequency domain method, a direct time-domain method such as might be programmed in C or Java, and a vectorized version such as would be common in Matlab or ...

8

The functions do not have a finite area, so they cannot be real distributions as your title claims they are. Let's change them a bit so they have area 1. f[x_] = (1/k) Exp[-x/k] UnitStep[x]; g[x_] = (1/p) Exp[-x/p] UnitStep[x]; Integrate[f[x], {x, -∞, ∞}] ConditionalExpression[1, Re[1/k] > 0] The convolution: Convolve[f[x], g[x], x, y] ...

7

While it's tempting to attribute the errors you're observing to floating-point errors due to zeros in the DFT of the window, this is actually not the case here: Window[width_, x_] := UnitStep[x + width/2] UnitStep[width/2 - x]; Test[x_] := UnitStep[x] UnitStep[count - x] Sin[1/10 x]^2; tempWindow = Table[Window[20, x], {x, -count/2, count/2}]; ...

6

Here is a way to solve this problem using the convolution theorem: l = Assuming[{\[Gamma] > 0 && \[Sigma] > 0 && \[Mu] > 0 && k \[Element] Reals}, FourierTransform[PDF[CauchyDistribution[\[Mu], \[Gamma]], x], x, k] ] $\frac{\left(\theta (-k) e^{2 \gamma k}+\theta (k)\right) e^{-k (\gamma -i \mu )}}{\sqrt{2 ... 6 I think you are trying to do the following: Clear[h]; h[n_] /; n >= 1 := -((I^-n (-1 + I^n)^2)/(n^2 Pi^2)) h[0] = 1/4; Clear[y]; y[m_] := DiscreteConvolve[h[n], DiscreteDelta[n], n, m] y[0] (* ==> 1/4 *) Here corrected your definition of y so it uses m as the function argument, but the important part is to understand why even with this ... 6 I'm not sure you are handling the top and left edges in the way you really want; they work with the second rather than first elements being treated as "middle" twice, with first elements not treated that way at all. Here is code that does not do that, hence gives different results than yours on top and left edges. It is around two orders of magnitude ... 5 The ideas mentioned in comments and the prior response seem like good ways to go about this. As for the brute force direct method, for a reliable result you can precompute one part symbolically and handle the rest numerically. ii[y_] = Integrate[ PDF[NormalDistribution[0, 8/1000], x - y]*4 *DiracDelta[1 - x], {x, 0, 11/10}]; firstTry[y_?NumberQ] := ... 5 Sparse multiplication with a circular matrix corresponds to a convolution; on a trivial example let us compare: matrix = SparseArray[{Band[{1, 1}] -> 2, Band[{1, 2}] -> 1, Band[{2, 1}] -> 1}, {15, 15}]; vec = SparseArray[5 -> x, 15]; matrix.vec // Normal (* ==> {0, 0, 0, x, 2 x, x, 0, 0, 0, 0, 0, 0, 0, 0, 0} *) versus a = SparseArray[5 ... 4 Using convolution theorem http://en.wikipedia.org/wiki/Convolution_theorem that Fourier of the convolution of two functions is the same as multiplications of their Fourier transforms then Use UnitStep to generate the time limited sin function to convolve with, like this Plot[UnitStep[t] UnitStep[Pi - t] Sin[t], {t, -3 Pi, 3 Pi}] and now apply the ... 4 I think the problem is that it could not do the symbolic integration. Meanwhile, you could always do the convolve directly, since it is just an integration, using NIntegrate. Unless you are looking for a closed form expression of conv(f,g). It would be more efficient to find the conv(f,g) expression, and then evaluate it for different$y$values than having ... 4 You have to go line by line to see your errors. You have to do this for each RGB channel separately. img = Import["ExampleData/lena.tif"]; data = ImageData /@ ColorSeparate[img]; {row, col} = Dimensions[data[[1]]]; H[a_, b_, T_] := Function[{x, y}, z = Pi*(a*x + b*y); T*N[Sinc[z]*Exp[-I*z]]] filter = Array[H[0.1, 0.1, 1], {row, col}]; filter // Abs // ... 4 As Nasser has observed numerical integration works. The following plots (f[x],f[y-x] f[x]f[y-x[ and convolution below) for varying y and thus gives insight into convolution). where convolution was obtained: g[y_] := NIntegrate[f[x] f[y - x], {x, -Infinity, Infinity}] 3 You can do the cross-correlation between two sequences using ListCorrelate a = {1, 2, 3, 4, 5}; b = {5, 4, 3, 2, 1}; ListCorrelate[a, b, 1] This does circular correlation, so you may want to look at the options to get the exact calculation you are looking for. 3 Try this: a.RotateLeft[b, 1] (*45*) For any other time delay$k$, use a.RotateLeft[b, k] instead. If the lists are extremely long and you are trying to compute the cross-correlation with a variety of delays, you can look into an FFT-based convolutive method. 3 Hard to do it analytically. Tried convolution theorem also. ForuierTransform had hard time with it as well as Integrate. So, here is a numerical solution. The support needed is really only from$0$to$2 \sqrt(2)$since your function exist over$0$to$\sqrt(2)\$ but I integrated it over little larger range for the plot to look better. Hence ...

3

You can write the definition down like that in Mathematica: r[0, t_] := Exp[-k t] Cos[ω0 t] r[n_, t_] := r[n, t] = Integrate[r[0, t - tt] r[n - 1, tt], {tt, 0, t}, Assumptions -> k > 0 && ω0 > 0] ... or is that what you considered to be brute force, and you're looking for a general solution of the recursion relation? ...

3

I'd suggest doing something like this Clear[x, y, c, d, exp] x[t_] := Exp[-2*t]*(HeavisideTheta[t + 2] - HeavisideTheta[t]) + (1 - t/2)*(HeavisideTheta[t] - HeavisideTheta[t - 2]) y[t_] := -2*DiracDelta[t + 2] + 2*(HeavisideTheta[t] - HeavisideTheta[t - 2]) c[t_] := (c[t] = Convolve[x[d], y[d], d, t]); kkkk = 4; arr = Array[c, kkkk*8 + 1, {-4, ...

3

Something like the following: dist = ProbabilityDistribution[ Convolve[PDF[NormalDistribution[4, 5], x], PDF[NormalDistribution[3, 1], x], x, t], {t, -Infinity, Infinity}] Expectation[z, z \[Distributed] dist] (* 7 *)

2

You can get some speed increase by writing your functions differently. It is almost always faster to act on whole lists in Mathematica than to explicitly loop over the individual elements. A couple of notes: Round is listable and does not have to be mapped over the array. The definition of tmp is just a dot product. The alterations to tmp in your If ...

Only top voted, non community-wiki answers of a minimum length are eligible