Overlapped stacking

[Guest]

This is for the purposes of doing photometry on a short period variable star.  Suppose I take 100 x 30 second images with say 10 second gap between the shots. I need to stack to get a reasonable SNR. I had this idea of doing a kind of overlapped stack.  Stacking 10 images at a time but each group of ten overlaps the previous set of ten. Like this -

Stack 1..10

Stack 5..14

Stack  10..19

Stack 15..24

Stack 20..29

And so on.

What I have done in the past is to stack in groups of 10. So 1-10, 11-20, 21-30 and so on.  In that case my light curve would have 10 points. But if I did the overlap stacking I'd have 20 points.   You could go as far as stacking in tens every image ie 1-10, 2-11, 3-12 and so on in which case I'd get 100 points on my light curve.

Is that a daft idea?

Cheers

Steve

 

[vlaiv]

Too much work as you can get the same result in easier way.

By doing that you are effectively interpolating between original data points. There is no more data to be gained - these intermediate stacking steps will have high correlation with other samples - so you might as well just interpolate existing 10 points with some interpolation algorithm to get nice curve.

check out this for example (hit example on the page and play around with data)>

https://tools.timodenk.com/cubic-spline-interpolation

[JRWASTRO]

Greetings Steve,

Re:  Is that a daft idea?

Definitely not ! What you are doing is a moving average. This is equivalent to low pass filtering where the noise bandwidth is determined by the bandwidth of this LPF. The more you "stack" the lower will be the bandwidth. But remember too much averaging will "smear" the results you will be looking for.

You are trying to obtain a graph of the brightness vs time so the sampling rate must be at least twice the highest rate of the brightness change - please see Nyquists sampling theorem.

You need to take into account the period of the changing brightness.   ie. the sampling frequency is 1/40 Hz.  If the change in brightness is very slow compared to 400s (ie. 10 x( 30s + 10s) ) then you are good to go!! ie take 30s samples and average 10 at a time as you have said. Your last strategy would be the way to go.

Jeremy.

 

[Guest]

Thank you both, rolling average was what I had in mind.  I had to cut short my observations last night because of cloud so I wasn't able to do a full cycle of about 98 minutes. The star was a cataclysmic variable which was probably not a good one to try.

Steve

 

[Xilman]

1 hour ago, vlaiv said:

Too much work as you can get the same result in easier way.

By doing that you are effectively interpolating between original data points. There is no more data to be gained - these intermediate stacking steps will have high correlation with other samples - so you might as well just interpolate existing 10 points with some interpolation algorithm to get nice curve.

check out this for example (hit example on the page and play around with data)>

https://tools.timodenk.com/cubic-spline-interpolation

Also look into Savitzky-Golay filtering. That method uses a least-squares polynomial fit to runs of data points. It works very well in practice in reducing noise without greatly broadening genuine features.

A simple variation can be used to find maxima and minima in the data, perhaps from an eclipsing binary.

Descriptions and code for S-G filtering are all over the net. Finding them is left as an exercise. Hint: Wikipedia is a good starting point.

[Guest]

Talking about Nyquist and filters. A long, long time ago I was involved in testing the frequency response of jet engines and associated fuel systems. Digitally sampling recordings of various engine parameters. One of the things we had to do was to use antialiasing filters before the sampling in order to avoid the fold-over effect.

 

[JRWASTRO]

Greetings Steve,

Implement the SG filter and presenting the results will need some software like MATLAB (very expensive) or OCTAVE (free open source software).

Create a list of your raw data (CSV) that can be read by OCTAVE and it will be straight forward to implement a SG filter and or a Moving Average (MA) filter.

For the SG filter you will need to calculate the polynomial coefficients at each point. Calculation of the coefficients involve some matrix manipulations like  multiplications and inverse of a matrix.  Octave is suited ideally to do these calculations.

Once you have obtained these polynomial coefficients a filtering algorithm may we written. This algorithm may  be used to do a MA filter (set all the coefficients = 1.0) or the SG filter. The SG filter is a general version of a MA filter and will be superior in preserving peaks.

Plotting / presentation of the results is presented easily  using OCTAVE. --  OCTAVE 6.2.0,  c.2021, is available for download.

Jeremy.

 

[vlaiv]

Moving average is rather "strange" low pass filter.

Fourier transform of box is Sinc function. Frequency response of moving average looks like this:

(colors show different size moving average)

More details here:

https://ptolemy.berkeley.edu/eecs20/week12/freqResponseRA.html

[JRWASTRO]

Nothing strange about this.  Interesting comparison though.

Re: We can do much better than that.

The author talks about a 16 moving average but does not tell us about the length of the impulse response of the LPF. Perhaps he/she needs to show the response of different length LPF's. 

 

[vlaiv]

9 hours ago, JRWASTRO said:

Nothing strange about this.  Interesting comparison though.

Re: We can do much better than that.

The author talks about a 16 moving average but does not tell us about the length of the impulse response of the LPF. Perhaps he/she needs to show the response of different length LPF's. 

 

I guess you are right.

I found it strange due to my expectations that low pass filter should be monotonically decreasing in frequency domain.

[Martin Meredith]

26 minutes ago, vlaiv said:

I guess you are right.

I found it strange due to my expectations that low pass filter should be monotonically decreasing in frequency domain.

When you think about it it makes sense because the abrupt edge of a rectangular window cannot ever be approximated as the sum of only low frequency components (deconstructing what the Fourier theorem). Try generating a sinusoid at some frequency f, then add in the first harmonic ie f + 2f, then add in also the next harmonic, etc and see what you get. Eventually, with high enough harmonics, you start to approximate a hard edge. (This assumes you add them all in zero phase, otherwise all bets are off but you won't get a hard edge either). Of course, the amplitude of the harmonic matters too and the above demonstration should be done with equal amplitudes to see the effect.

Here's a nice discussion of the spectral effect of various window functions: https://en.wikipedia.org/wiki/Window_function#Windowing

I forget where I heard it, but it was once remarked that practically any function that rises to a peak and then falls again has someone's name on it as a window function...

Martin

Edited March 26, 2021 by Martin Meredith

[vlaiv]

6 minutes ago, Martin Meredith said:

When you think about it it makes sense because the abrupt edge of a rectangular window cannot ever be approximated as the sum of only low frequency components (deconstructing what the Fourier theorem). Try generating a sinusoid at some frequency f, then add in the first harmonic ie f + 2f, then add in also the next harmonic, etc and see what you get. Eventually, with high enough harmonics, you start to approximate a hard edge. (This assumes you add them all in zero phase, otherwise all bets are off but you won't get a hard edge either).

Here's a nice discussion of the spectral effect of various window functions: https://en.wikipedia.org/wiki/Window_function#Windowing

I forget where I heard it, but it was once remarked that practically any function that rises to a peak and then falls again has someone's name on it as a window function...

Martin

Oh it makes perfect sense if you mean the shape of FT of box function - that part I'm fine with. It's just the way FT calculates and yes - box and sinc are know dual of FT.

It is just usability of it as low pass filter that I was referring to. For some strange reason I expected a good low pass filter to be monotonically decreasing in frequency domain - like say Gaussian blur filter, but in reality, why should it. I guess that only requirement for filter to be called low pass is to attenuate higher frequencies more than lower, right?

[Xilman]

12 hours ago, JRWASTRO said:

Greetings Steve,

Implement the SG filter and presenting the results will need some software like MATLAB (very expensive) or OCTAVE (free open source software).

Create a list of your raw data (CSV) that can be read by OCTAVE and it will be straight forward to implement a SG filter and or a Moving Average (MA) filter.

Octave is a good program but I dispute the necessary part. Anyone who followed my advice and started investigations with the Wikipedia article would have seen that coefficients up to degree 25 were published as long ago as 1965.

A link in the wikipedia leads to https://sites.google.com/site/chandraacads/resources/sg-filter/db

where you will find tables of coefficients for 1 through 4 dimensional data and for a wide range of polynomial degrees and window sizes.

Once you have the coefficients it is relatively easy to implement the filter in the programming language of your choice, Excel-compatible spreadsheets if you wish though, to be honest in high dimensions it will run rather slowly unless you use FFT convolution.

That said, if the tabulated coefficients are not suitable for your needs, you will need to compute them. Free software exists to do this for you. Another link from the Wikipedia page leads to C source https://zenodo.org/record/1288901

Have fun!

 

[JRWASTRO]

Very good.

Please tell me how many polynomial coefficient sets do I need to use here:

My noisy data set has c. 22000 points  of a complex waveform. The data is, of necessity, complex it is at 0Hz.  LPF here produces a BPF.

Jeremy.

[Xilman]

8 hours ago, JRWASTRO said:

Very good.

Please tell me how many polynomial coefficient sets do I need to use here:

My noisy data set has c. 22000 points  of a complex waveform. The data is, of necessity, complex it is at 0Hz.  LPF here produces a BPF.

Jeremy.

I don't know is the simple answer. The reason is that I have only ever applied S-G filtering to purely real data. Complex data is outside my experience. Is its Fourier transform pure real?

The original question was for a real time series and my answer was given in that context.

However, a different and interesting question which I may investigate further.

[JRWASTRO]

Not to worry, as we say.

I only wish to point out that the filtering algorithm should calculation of the polynomial coefficients at each step. In the example I used if you use a polynomial of order 11 and padded the data file with 5 points on the front and back one will need to compute 2000 polynomial sets.

About the DFT  - it is evaluated using the Fast Fourier Transform - but that is another issue!!

The FT of a signal is complex - it has a magnitude and a phase. For real signals the FT will exhibit Hermitian Symmetry (a very important property). When you multiply the X(f) by its complex conjugate X*(f) you will get a real value representing a power spectrum etc. etc. etc.

FYI I can take two signals x1(t)  => X1(f) and x2(t) => X2(f) and obtain : S12(f) = X1(f) x X2*(f) which is the cross power spectral density and is a complex function which tells us about the correlation between the two functions - especially the phase of the PSD.

Jeremy.

 

 

[Xilman]

On 28/03/2021 at 00:07, JRWASTRO said:

Not to worry, as we say.

I only wish to point out that the filtering algorithm should calculation of the polynomial coefficients at each step. In the example I used if you use a polynomial of order 11 and padded the data file with 5 points on the front and back one will need to compute 2000 polynomial sets.

About the DFT  - it is evaluated using the Fast Fourier Transform - but that is another issue!!

The FT of a signal is complex - it has a magnitude and a phase. For real signals the FT will exhibit Hermitian Symmetry (a very important property). When you multiply the X(f) by its complex conjugate X*(f) you will get a real value representing a power spectrum etc. etc. etc.

FYI I can take two signals x1(t)  => X1(f) and x2(t) => X2(f) and obtain : S12(f) = X1(f) x X2*(f) which is the cross power spectral density and is a complex function which tells us about the correlation between the two functions - especially the phase of the PSD.

Jeremy.

 

 

For your first point, the beauty of the S-G algorithm is that the polynomial sets are computed iteratively and represented as a convolution kernel. For those readers who have not yet learned the theory behind it, I urge you to do so. It really is a beautiful algorithm and the original paper is deservedly one of the most quoted in the scientific literature.

I know about how to compute Fourier transforms, their properties, and how to use them for various convolution, correlation, spectral analysis and so on. I kept my presentation simple and to the subject of smoothing of real valued time series data. Anyone who wishes to learn more can easily find vast amounts of information on the web, ranging from hand-waving to hard-core mathematics and everything in between, including software implementations.

In your particular case, wherein lies the significant noise? x1(t) or x2(t) or both?  Have you characterized the noise (e.g, pink, white or blue)? Have you tried smoothing one or both time series? Sorry I can't help much but that is because I have never needed to address the problem of interest. Presumably you have done a literature search.

Edited March 31, 2021 by Xilman