I'm working on a "problem" that got me thinking.
As the purpose of attempting the "problem" was education that's good.
I have 96 data points.
I desire:
1. approximation of 97th (but NOT 98th or greater).
2. interpolate 28 intervening data points, including between points
96 and 97.
I know a priori that:
1. the function and ALL derivatives are continuous.
2. the function is almost periodic.
a. 1 period ends ~3/4 ow way between points 96 and 97.
b. The deviations from periodic behavior are small and in
themselves periodic but with periods 30 to 300 times longer
than my total sample time.
Scilab's smooth() function appears to handle the interpolations from
points 1 thru 96.
Can DSP techniques or point of view:
1. get me a better interpolation?
2. get me point 97?
[ The 96 data points are satellite positions in an Earth Centered Earth
Fixed (ECEF) frame of reference from Rinex formated data files running
nominally from midnight to midnight. Unfortunately, the closing
"midnight" is in the next day's data file (which will not always be
accessible). As orbits are well defined, I could convert from ECEF to
inertial frames, interpolate and extrapolate, and convert back to ECEF.
That has it's own set of *MASSIVE* problems for me.]
Comments?
Interpolation and extrapolation of almost periodic data
Started by ●April 12, 2008
Reply by ●April 12, 20082008-04-12
On Apr 12, 9:33�am, Richard Owlett <rowl...@atlascomm.net> wrote:> I'm working on a "problem" that got me thinking. > As the purpose of attempting the "problem" was education that's good. > > I have 96 data points. > I desire: > � �1. approximation of 97th (but NOT 98th or greater). > � �2. interpolate 28 intervening data points, including between points > � � � 96 and 97. > > I know a priori that: > � �1. the function and ALL derivatives are continuous. > � �2. the function is almost periodic. > � � � a. 1 period ends ~3/4 ow way between points 96 and 97. > � � � b. The deviations from periodic behavior are small and in > � � � � �themselves periodic but with periods 30 to 300 times longer > � � � � �than my total sample time. > > Scilab's smooth() function appears to handle the interpolations from > points 1 thru 96. > > Can DSP techniques or point of view: > � �1. get me a better interpolation? > � �2. get me point 97? > > [ The 96 data points are satellite positions in an Earth Centered Earth > Fixed (ECEF) frame of reference from Rinex formated data files running > nominally from midnight to midnight. Unfortunately, the closing > "midnight" is in the next day's data file (which will not always be > accessible). As orbits are well defined, I could convert from ECEF to > inertial frames, interpolate and extrapolate, and convert back to ECEF. > That has it's own set of *MASSIVE* problems for me.] > > Comments?this is the kinda thing i think about regarding synthesis/analysis of harrmonic (or "quasi-periodic") musical tones. but i usually think i need about 2 adjacent periods to be able to perform some sorta pitch detection. if you know the precise period (to less that a sample precision), you ougta be able to interpolate, but the edge effects are a problem. r b-j
Reply by ●April 12, 20082008-04-12
Richard Owlett wrote:> > I'm working on a "problem" that got me thinking. > As the purpose of attempting the "problem" was education that's good. > > I have 96 data points. > I desire: > 1. approximation of 97th (but NOT 98th or greater). > 2. interpolate 28 intervening data points, including between points > 96 and 97.I don't understand what number 2 means. Are you increasing the data from 96 evenly spaced samples to 125 evenly spaced samples? Or are you reducing it to 27 samples? Anyway if you are using an interpolation method to create a point between point 95 and 96 that probably means that a 97th point was already invented to accomplish that. That is, if you are doing anything more than linear interpolation then to find points in the interval between 95 and 96 you need points outside that interval. You say you know the data is periodic. Wouldn't that mean that the missing point 97 is actually the same as point 1 ? -jim> > I know a priori that: > 1. the function and ALL derivatives are continuous. > 2. the function is almost periodic. > a. 1 period ends ~3/4 ow way between points 96 and 97. > b. The deviations from periodic behavior are small and in > themselves periodic but with periods 30 to 300 times longer > than my total sample time. > > Scilab's smooth() function appears to handle the interpolations from > points 1 thru 96. > > Can DSP techniques or point of view: > 1. get me a better interpolation? > 2. get me point 97? > > [ The 96 data points are satellite positions in an Earth Centered Earth > Fixed (ECEF) frame of reference from Rinex formated data files running > nominally from midnight to midnight. Unfortunately, the closing > "midnight" is in the next day's data file (which will not always be > accessible). As orbits are well defined, I could convert from ECEF to > inertial frames, interpolate and extrapolate, and convert back to ECEF. > That has it's own set of *MASSIVE* problems for me.] > > Comments?----== Posted via Pronews.Com - Unlimited-Unrestricted-Secure Usenet News==---- http://www.pronews.com The #1 Newsgroup Service in the World! >100,000 Newsgroups ---= - Total Privacy via Encryption =---
Reply by ●April 12, 20082008-04-12
On Apr 12, 6:43�pm, jim <"sjedgingN0sp"@m...@mwt.net> wrote:> > � � � � You say you know the data is periodic. Wouldn't that mean that the > missing point 97 is actually the same as point 1 ?i thought what happened is that he knows from some other source, what is the period. and that it has a fractional part, in terms of samples periods. but it is not sampled synchronously. the period might be between 95 and 96 samples long. that's what i thought the problem is. if it is, isn't that a matter of setting up some Fourier series with the known period and with no more than 47 harmonics + DC to fit those 96 points? other than that guess, i'm clueless to what the problem is. r b-j
Reply by ●April 12, 20082008-04-12
On Apr 12, 6:33 am, Richard Owlett <rowl...@atlascomm.net> wrote:> I'm working on a "problem" that got me thinking. > As the purpose of attempting the "problem" was education that's good. > > I have 96 data points. > I desire: > 1. approximation of 97th (but NOT 98th or greater). > 2. interpolate 28 intervening data points, including between points > 96 and 97. > > I know a priori that: > 1. the function and ALL derivatives are continuous. > 2. the function is almost periodic. > a. 1 period ends ~3/4 ow way between points 96 and 97. > b. The deviations from periodic behavior are small and in > themselves periodic but with periods 30 to 300 times longer > than my total sample time. > > Scilab's smooth() function appears to handle the interpolations from > points 1 thru 96. > > Can DSP techniques or point of view: > 1. get me a better interpolation? > 2. get me point 97? > > [ The 96 data points are satellite positions in an Earth Centered Earth > Fixed (ECEF) frame of reference from Rinex formated data files running > nominally from midnight to midnight. Unfortunately, the closing > "midnight" is in the next day's data file (which will not always be > accessible). As orbits are well defined, I could convert from ECEF to > inertial frames, interpolate and extrapolate, and convert back to ECEF. > That has it's own set of *MASSIVE* problems for me.] > > Comments?It's not clear just what you want, but we never let that stop us on comp.dsp. Assuming uniformly spaced samples, there are DFT based methods for decomposing signals into component parts. For signals that are accurately decomposable into a sum of sinusoids, there are DFT based decompositions. With the frequency phase and amplitude of all components, the signal can be calculated for all time subject to limitations of SNR, stationarity and model match. If the components are adequately separated, the process works as follows: 1) Perform DFT of the data set 2) Determine the coarse location of the strongest peak by peak-picking the magnitude of the DFT 3) Estimate the fine resolution frequency, amplitude and phase of the peak signal from the DFT coefficients of the bin with greatest magnitude and the coefficients of the adjacent bins. 4) Use the frequency, amplitude and phase of the peak to remove its contribution to the signal. In the frequency domain this includes all leakage effects, in the time domain it is the generation and subtraction of the estimated signal component. 5) When the residual energy is small continue, else go back to 1) or 2) as required by the domain used for component removal. 6) Use the frequencies, amplitudes and phases of the components to estimate the value of the components at any desired time and sum the components for the estimate. A discussion of the process is found in: M. D. Macleod, \Fast nearly ML estimation of the parameters of real or complex single tones or resolved multiple tones," IEEE Transactions of Signal Processing, vol. 46, pp. 141{148, Jan. 1998. An example of the application of the process to measuring superimposed tones in an intercept receiver you can download: ADA426470.pdf proxy Url: http://handle.dtic.mil/100.2/ADA426470 Title: Enhancing the Instantaneous Dynamic Range of Electronic Warfare Receivers Using Statistical Signal Processing Personal Author(s): Smith, Bryan E. Report Date: MAR 2004 Media Count: 121 Pages(s) This thesis includes a scheme for removing leakage errors in the frequency domain. If you want to do more computation and claim a more exact maximum likelihood approach download: ADA310879.pdf proxy Url: http://handle.dtic.mil/100.2/ADA310879 Title: Parameter Estimation for Superimposed Weighted Exponentials. Fields and Groups : 120300 - STATISTICS AND PROBABILITY 120400 - OPERATIONS RESEARCH 201400 - RADIOFREQUENCY WAVE PROPAGATION Corporate Author: AIR FORCE INST OF TECH WRIGHT-PATTERSON AFB OH SCHOOL OF ENGINEERING Personal Author(s): Ingham, Edwards A. Report Date: JUL 1996 Media Count: 212 Pages(s) If constant frequency sinusoids are not a good match for your signal (if the residuals in step 5 don't get small it's a hint) you can use other signal models such as linear FMs or polynomial phase signals at greater computational cost. Accuracy can be improved by reestimating the parameters of the larger peaks with the other peaks removed, it's a just a lot more computation. Some signals may require it. There is no need to limit the decomposition to DFT based methods, but if they work they may be computationally advantageous. If the samples are not uniformly spaced in time the Lomb-Scargle periodogram algorithms can be used to generate Fourier coefficients that are at equi-spaced frequencies. Dale B. Dalrymple http://dbdimages.com
Reply by ●April 13, 20082008-04-13
On Apr 12, 6:33 am, Richard Owlett <rowl...@atlascomm.net> wrote:> I'm working on a "problem" that got me thinking. > As the purpose of attempting the "problem" was education that's good. > > I have 96 data points. > I desire: > 1. approximation of 97th (but NOT 98th or greater). > 2. interpolate 28 intervening data points, including between points > 96 and 97. > > I know a priori that: > 1. the function and ALL derivatives are continuous.Do you have a model of the system producing the signal or function? Was the function or signal band-limited before the sampling? If so, what was the bandwidth? Do you know enough about the function to estimate what portion of the sample content is noise vs. signal?> 2. the function is almost periodic. > a. 1 period ends ~3/4 ow way between points 96 and 97.How many total periods do you have represented in the 96 point vector?> b. The deviations from periodic behavior are small and in > themselves periodic but with periods 30 to 300 times longer > than my total sample time. > > Scilab's smooth() function appears to handle the interpolations from > points 1 thru 96. > > Can DSP techniques or point of view: > 1. get me a better interpolation? > 2. get me point 97?If the noise level of the signal and it's derivatives are zero, then you might be able to use some number of finite differences (or other polynomial method) to interpolate forward. If the signal is noisy but all in narrow band(s), you might be able to filter it, or model it with an "adapted" all pole filter. If you have a parametric model of the system producing the signal, then you might be able to do some sort of linear (least squares?) or non-linear (genetic evolution?) best fit of the parameters to the samples, and then run the model forward for one more time step. If the signal consists of only one fundamental period, but has some rich harmonic content, then you could try a combination of changing the fft size and resampling the signal until the fundamental and every harmonic all fall at exact bin centers. (If the harmonics are known to be exact, then this may converge faster than dbd's suggestion of trying to remove each frequency peak one-at-a-time, since you will be using the information from several peaks to estimate only two parameters.) IMHO. YMMV. -- rhn A.T nicholson d.0.t C-o-M
Reply by ●April 13, 20082008-04-13
robert bristow-johnson wrote:> On Apr 12, 6:43 pm, jim <"sjedgingN0sp"@m...@mwt.net> wrote: > >> You say you know the data is periodic. Wouldn't that mean that the >>missing point 97 is actually the same as point 1 ? > > > i thought what happened is that he knows from some other source, what > is the period. and that it has a fractional part, in terms of samples > periods. but it is not sampled synchronously. the period might be > between 95 and 96 samples long.Yes to knowing period from other sources. Points 1 -> 97 cover exactly 24 hours. The period is 23 hours 56 minutes and some seconds. Point 97 is "available" in another file, but that file may not be accessible. There are workarounds for my specific case. *BUT* the problem got me thinking about the general case. THEREFORE I posed the question in my original post with NO explicit reference to the specific case. The *PURPOSE* was to educate me on how to _THINK_ about DSP problems.> > that's what i thought the problem is. if it is, isn't that a matter > of setting up some Fourier series with the known period and with no > more than 47 harmonics + DC to fit those 96 points? > > other than that guess, i'm clueless to what the problem is. > > r b-j
Reply by ●April 13, 20082008-04-13
The bon-n-n-g-g-g you heard this morning was from your 100 lb sledge shakin out some cobwebs ;) In the mid 60's I maintained what might be considered a crude analog computer that used essentially your approach. dbd wrote:> On Apr 12, 6:33 am, Richard Owlett <rowl...@atlascomm.net> wrote: > >>I'm working on a "problem" that got me thinking. >>As the purpose of attempting the "problem" was education that's good. >> >>I have 96 data points. >>I desire: >> 1. approximation of 97th (but NOT 98th or greater). >> 2. interpolate 28 intervening data points, including between points >> 96 and 97. >> >>I know a priori that: >> 1. the function and ALL derivatives are continuous. >> 2. the function is almost periodic. >> a. 1 period ends ~3/4 ow way between points 96 and 97. >> b. The deviations from periodic behavior are small and in >> themselves periodic but with periods 30 to 300 times longer >> than my total sample time. >> >>Scilab's smooth() function appears to handle the interpolations from >>points 1 thru 96. >> >>Can DSP techniques or point of view: >> 1. get me a better interpolation? >> 2. get me point 97? >> >>[ The 96 data points are satellite positions in an Earth Centered Earth >>Fixed (ECEF) frame of reference from Rinex formated data files running >>nominally from midnight to midnight. Unfortunately, the closing >>"midnight" is in the next day's data file (which will not always be >>accessible). As orbits are well defined, I could convert from ECEF to >>inertial frames, interpolate and extrapolate, and convert back to ECEF. >>That has it's own set of *MASSIVE* problems for me.] >> >>Comments? > > > It's not clear just what you want, but we never let that stop us on > comp.dsp.So consider the OP a unknown N-port device. So instead of hitting it with a pure impulse to get total transfer function, you hit it with a finite set of ideas.> > Assuming uniformly spaced samples,They are> there are DFT based methods for > decomposing signals into component parts. For signals that are > accurately decomposable into a sum of sinusoids, there are DFT based > decompositions. With the frequency phase and amplitude of all > components, the signal can be calculated for all time subject to > limitations of SNR, stationarity and model match.The system has 4 independent signal sources. Two of which are comparable in magnitude of their fundamentals. Their frequencies differ by ~.2% . The other are 30-300+ times slower and much weaker. But ALL 4 fundamental frequencies are known very accurately.> > If the components are adequately separated, the process works as > follows: > > 1) Perform DFT of the data set > > 2) Determine the coarse location of the strongest peak by peak-picking > the magnitude of the DFT > > 3) Estimate the fine resolution frequency, amplitude and phase of the > peak signal from the DFT coefficients of the bin with greatest > magnitude and the coefficients of the adjacent bins. > > 4) Use the frequency, amplitude and phase of the peak to remove its > contribution to the signal. In the frequency domain this includes all > leakage effects, in the time domain it is the generation and > subtraction of the estimated signal component. > > 5) When the residual energy is small continue, else go back to 1) or > 2) as required by the domain used for component removal. > > 6) Use the frequencies, amplitudes and phases of the components to > estimate the value of the components at any desired time and sum the > components for the estimate. > [snip references and discussion of non-uniform sampling and varyingsignal frequencies] Got to get dressed for church. More later.
Reply by ●April 13, 20082008-04-13
> Yes to knowing period from other sources. > Points 1 -> 97 cover exactly 24 hours. > The period is 23 hours 56 minutes and some seconds.Aha a sidereal period. Richard, there are quite a few sources that will let you find the osculating (not oscillating!) coefficients. A good book on orbital mechanics will likely prove useful. Given your simple case of data and limited interval for extrapolation, you can try a simple LPC method. It should be good enough. Clay
Reply by ●April 13, 20082008-04-13
clay@claysturner.com wrote:>>Yes to knowing period from other sources. >>Points 1 -> 97 cover exactly 24 hours. >>The period is 23 hours 56 minutes and some seconds. > > > Aha a sidereal period. Richard, there are quite a few sources that > will let you find the osculating (not oscillating!) coefficients. A > good book on orbital mechanics will likely prove useful.I've a limited budget (read as ~null ;) What would be good keywords for Google searching?> > Given your simple case of data and limited interval for extrapolation, > you can try a simple LPC method. It should be good enough. >Yep. And Dale used a large enough hammer to wake me up to already "knowing" the answer ;)> Clay >
