2
$\begingroup$

I am trying to rationalize a figure given in the Numerical Recipes in C in the section of Fourier based convolution and deconvolution.

The authors show the example of a continuous convolution with a zero centered response function r(t). Convoluting the continuous signal blurs out the edges and there is no shift.

What is somewhat confusing is that they suggestion for the discrete convolution analogue of the continuous convolution, we should wrap the negative time section of the response function to the extreme right end of the array. What is the mathematical justification of doing N/2 circular shift, where N is even number of samples. The authors do not mention as N/2 circular shift, but it seems this is what they are doing.

These are their actual figures.

enter image description here

It seems the authors are suggesting circular shifting the zero centered response by N/2 (N is even) to get an equivalent convolution to continuous convolution. Why can't we have discrete version of the response function as shown by the squared markers which I drew in MATLAB? I do not get equivalent convolution results if we convolute a square wave with the sampled response function on the right and the left as suggested in Numerical Recipes.

enter image description here

$\endgroup$
4
  • $\begingroup$ are you aware of the convolution property of the Fourier Transform? And of the Discrete Fourier Transform? $\endgroup$ Jan 22 at 8:11
  • 1
    $\begingroup$ Yes, are you referring to the multiplication property of the continuous function by continuous FT and the circular nature of the DFT convolution? $\endgroup$
    – M. Farooq
    Jan 22 at 14:07
  • $\begingroup$ My main question is how the continuous convolution of Fig 13.1.1 is equivalent to discrete convolution of Fig 13.1.2 ? $\endgroup$
    – M. Farooq
    Jan 22 at 14:10
  • $\begingroup$ you just answered that yourself, I think? $\endgroup$ Jan 22 at 14:18
2
+50
$\begingroup$

It's about zero-phase filtering.

First, the continuous-time response function $r(t)$ is an even-symmetric, noncausal filter so that its Fourier transform will be a zero-phase function; i.e., $R(\omega)$ is a real-valued function of frequency $\omega$.

Note that there's almost a mistake in the figure 13.1.1; the output should begin before $t=0$, even for a zero-phase filter. Because this is a non-causal filter, it will anticipate future inputs before they happen, and put out the response before $t=0$. I think, however, they wanted to justify this condition, by shifting the input $s(t)$ sufficiently to the right (or assuming that $s(t)=0$ until $t=t_1$) so that the output response before $t=0$ will be zero anyway, and thus eventually we are interested at the output of the filter beginning from $t=0$.

Now, sample the signals $s(t)$ and $r(t)$ (assuming they are bandlimited) and obtain a sampled version of the continuous-time output by performing a discrete-time linear convolution, or equivalently use a DFT based efficient implementation of linear convolution by circular convolution of proper length.

If you choose the first approach, and perform direct discrete-time linear convolution using the sequences $s[n]$ and $r[n]$, then you won't need to shift anything. But remember that a practical implementation will still need some shifts to represent the noncausal impulse response $r[n]$ in the proper index frame.

But if you choose the second approach; using DFT to implement a linear filtering, you have two choices for the filter response $r[n]$ to be registered; (because in the DFT domain, there's no negative time-index $n$, and you are restricted to the range $0 \leq n < N$) first you can linearly shift $r[n]$ sufficiently right to make it a causal filter $r_c[n]$, and follow the procedure as usual; with this choice, the effective discrete-time filter loses its zero-phase (and zero group delay) property and your effective output begins at the group delay position index.

If, instead, you wrap the negative-index portion of $r[n]$ to its far right edge, in the DFT frame of $N$ samples, where $N$ is chosen properly to aviod aliasing in the time-domain to ensure that computed circular convolution will match the desired linear one, then you will maintain the zero-phase property of $r[n]$ while performing the DFT based convolution; your first output at $n=0$ will be what you would get at $t=0$ in the continuous-time case.

$\endgroup$
8
  • $\begingroup$ Thanks. What is the mathematical process called is DSP when we wrap around the negative portion of the continuous version of the same function to maintain zero phase property? I could not find this discussed anywhere except this Numerical Recipe. $\endgroup$
    – M. Farooq
    Jan 29 at 1:05
  • $\begingroup$ The second point which is bothering me is that suppose we have a sequence, S= [0,0,0,1,1,0,0,0], if I linearly convolute it with a[n]= [0, 0, 0, 0.018, 0.106, 0.005, 0, 0] vs. b[n]=[0.106, 0.005, 0, 0, 0, 0,0, 0.018], I get very different results. How come a[n] and b[n] are equivalent? $\endgroup$
    – M. Farooq
    Jan 29 at 1:12
  • $\begingroup$ @M.Farooq I do not know a name given specifically to this operation. However, technically what you do is a circular left shift of the linearly right shifted sequence. You first linearly right shift the noncausal r[n] to make it causal, and then circularly right shift (by about half of its length) to make it effectively zero-phase (non-causal) again. It does not have a name. Another interpretation is that you naturally recast the noncausal r[n] into a DFT frame of [ 0 , N-1] , by interpreting r[n] as a periodic sequence with a period N. $\endgroup$
    – Fat32
    Jan 29 at 1:16
  • $\begingroup$ Thanks, do you remember a book which discusses these operations? Actually, I am interested in deconvolution of spectroscopic signals. What I am actually doing is that I have an experimental spectrum of N points. Then I use the re-casted version of a zero centered Lorentzian to prevent any position shift after dividing the DFT of a signal S, and that of $recasted$ Lorentzian. I found this by trying several ways so I was trying to rationalize as to why the recasted version does not change the position of peaks. $\endgroup$
    – M. Farooq
    Jan 29 at 1:23
  • 1
    $\begingroup$ these operations are distributed there and there but I do suggest DFT chapter of any standard DSP textbook, such as from Proakis, Oppenheim, Orfanidis (free), etc. $\endgroup$
    – Fat32
    Jan 29 at 1:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.