Forums

Discrete Deconvolution

Started by Unknown July 21, 2006
Hello, all.

I have a dilemma.  I'm analyzing some time signals for periodicities,
and I was wondering about a seemingly basic method of doing it.  I hope
you'll forgive my ignorance, as my specialty is physics and not signal
processing :P

So, we have a time-series that is the product of a source signal and a
window function :

c(t_i) = w(t_i)s(t_i)

where t_i represents the i_th time bin, c is the measured signal, w is
a window function, and s is the actual source signal.  I know the FFT
of c(t_i) = c(\omega_i) , and the FFT of the window function,
w(\omega_i).  However, I would really like to know the FFT of s,
s(\omega_i).

I hear about deconvolution from image analysis in astronomy and other
places, but am unfamiliar with it in time-series analysis.  I know it's
an ill-posed mathematical problem, but is there a way to extract
information about s(\omega_i) by using some sort of deconvolution
method?  Thanks.

Robert Harris

harrisrj@mit.edu wrote:

> Hello, all. > > I have a dilemma. I'm analyzing some time signals for periodicities, > and I was wondering about a seemingly basic method of doing it. I hope > you'll forgive my ignorance, as my specialty is physics and not signal > processing :P > > So, we have a time-series that is the product of a source signal and a > window function : > > c(t_i) = w(t_i)s(t_i) > > where t_i represents the i_th time bin, c is the measured signal, w is > a window function, and s is the actual source signal. I know the FFT > of c(t_i) = c(\omega_i) , and the FFT of the window function, > w(\omega_i). However, I would really like to know the FFT of s, > s(\omega_i). > > I hear about deconvolution from image analysis in astronomy and other > places, but am unfamiliar with it in time-series analysis. I know it's > an ill-posed mathematical problem, but is there a way to extract > information about s(\omega_i) by using some sort of deconvolution > method? Thanks. >
Do you no longer have the original time series? If you _have_ s(t_i) just take it's FFT. If you don't, then take the inverse FFT of c(t_i) and divide it by the inverse FFT of w(t_i) point by point, taking the fact that w(t_i) probably has some points equal to zero. But this is a bad idea for several reasons. In all but some very rare cases the FFT of a time series is an approximation of some other Fourier transform that you'd like to take but can't. The FFT, by itself, is just a really fast way of finding the Fourier series of a periodic waveform in discrete time, who's period is equal to the number of points you're feeding into the FFT. That's all. The FFT can be used as an approximation of the Fourier transform of a signal in sampled time with an infinite number of points, which is in turn an approximation of a signal in continuous time that extends into infinity. Like all approximations, however, the FFT has its problems when it's used in this way. The biggest problem is that while _you_ think you're taking the Fourier transform of a signal that happens to be of finite extent, the FFT thinks that it's taking the transform of a signal that repeats, that is that s(t_N) = s(t_0) for an N-point FFT. Usually this is a whopping big discontinuity, which makes the algorithm report lots of really high-frequency content which simply isn't there in the real signal. This big problem is solved by windowing, as a way to get the algorithm to cough up an approximate Fourier transform that has some hope of being close. I can almost guarantee you that if you take away the windowing you'll take away any hope of having an accurate representation of your original process. But you can find this out yourself: Pick an N, and generate a cosine wave that has a frequency of 10.5/N -- that is, it starts at 1 at t_0, and goes to nearly -1 at t_{N-1}. Now take a FFT of it -- does it look like a clean spike, balanced between bins 10 and 11? Now window it, and do another FFT -- any better? Perhaps if you tell us what you want to learn by all this we can help you find the right algorithm. -- Tim Wescott Wescott Design Services http://www.wescottdesign.com Posting from Google? See http://cfaj.freeshell.org/google/ "Applied Control Theory for Embedded Systems" came out in April. See details at http://www.wescottdesign.com/actfes/actfes.html
Tim Wescott wrote:

   ...

> The biggest problem is that while _you_ think you're taking the Fourier > transform of a signal that happens to be of finite extent, the FFT > thinks that it's taking the transform of a signal that repeats, that is > that s(t_N) = s(t_0) for an N-point FFT. Usually this is a whopping big > discontinuity, which makes the algorithm report lots of really > high-frequency content which simply isn't there in the real signal. > > This big problem is solved by windowing, as a way to get the algorithm > to cough up an approximate Fourier transform that has some hope of being > close. I can almost guarantee you that if you take away the windowing > you'll take away any hope of having an accurate representation of your > original process. > > But you can find this out yourself: Pick an N, and generate a cosine > wave that has a frequency of 10.5/N -- that is, it starts at 1 at t_0, > and goes to nearly -1 at t_{N-1}. Now take a FFT of it -- does it look > like a clean spike, balanced between bins 10 and 11? Now window it, and > do another FFT -- any better? > > Perhaps if you tell us what you want to learn by all this we can help > you find the right algorithm.
Tim, you too should teach. I guess writing books comes close. Jerry -- Engineering is the art of making what you want from things you can get. �����������������������������������������������������������������������
Jerry Avins wrote:

> Tim Wescott wrote: > > ... > >> The biggest problem is that while _you_ think you're taking the >> Fourier transform of a signal that happens to be of finite extent, the >> FFT thinks that it's taking the transform of a signal that repeats, >> that is that s(t_N) = s(t_0) for an N-point FFT. Usually this is a >> whopping big discontinuity, which makes the algorithm report lots of >> really high-frequency content which simply isn't there in the real >> signal. >> >> This big problem is solved by windowing, as a way to get the algorithm >> to cough up an approximate Fourier transform that has some hope of >> being close. I can almost guarantee you that if you take away the >> windowing you'll take away any hope of having an accurate >> representation of your original process. >> >> But you can find this out yourself: Pick an N, and generate a cosine >> wave that has a frequency of 10.5/N -- that is, it starts at 1 at >> t_0, and goes to nearly -1 at t_{N-1}. Now take a FFT of it -- does >> it look like a clean spike, balanced between bins 10 and 11? Now >> window it, and do another FFT -- any better? >> >> Perhaps if you tell us what you want to learn by all this we can help >> you find the right algorithm. > > > Tim, you too should teach. I guess writing books comes close. > > Jerry
That's on the agenda: look at the 'seminars' page of my web site for what I want to be doing next year. -- Tim Wescott Wescott Design Services http://www.wescottdesign.com Posting from Google? See http://cfaj.freeshell.org/google/ "Applied Control Theory for Embedded Systems" came out in April. See details at http://www.wescottdesign.com/actfes/actfes.html
<harrisrj@mit.edu> ??? 
??????:1153513373.755504.77180@m79g2000cwm.googlegroups.com...
> Hello, all. > > I have a dilemma. I'm analyzing some time signals for periodicities, > and I was wondering about a seemingly basic method of doing it. I hope > you'll forgive my ignorance, as my specialty is physics and not signal > processing :P > > So, we have a time-series that is the product of a source signal and a > window function : > > c(t_i) = w(t_i)s(t_i) > > where t_i represents the i_th time bin, c is the measured signal, w is > a window function, and s is the actual source signal. I know the FFT > of c(t_i) = c(\omega_i) , and the FFT of the window function, > w(\omega_i). However, I would really like to know the FFT of s, > s(\omega_i). > > I hear about deconvolution from image analysis in astronomy and other > places, but am unfamiliar with it in time-series analysis. I know it's > an ill-posed mathematical problem, but is there a way to extract > information about s(\omega_i) by using some sort of deconvolution > method? Thanks. > > Robert Harris
--------------------------------------------------------------------- Hi Convolution is some kind of filter. When you filter something you loose the data that stay outside the filter. The filter can be 100 db loose of signal and even more. In this case you have to "amplified back" 100 db of signal and for any practical reasons it is impossible. Many convolutiones can not go back and therefor deconvolution is an ill subject. The windows generally destroy the front and the back of the data and you have to make IFFT and then mutiply by the 1/window in order to have back the signal. If the window is low the 1/window is very high . I think that in many cases you can't do deconvolution at all. Best Regards Moshe
>