Guest post by Richard Saumarez
Correlations are often used to relate changes in climate variables in order to establish a potential relationship between them. These variables may have been processed in some way before they are correlated, possibly by filtering them. There has been some debate, which may have shed more heat than light, about the validity of combining these processes and whether they interact to make the conclusions drawn from them invalid. The object of this post is to explain the processes of correlating and filtering signals, the relationship between them and show that the results are predictable provided one takes care.
The importance of the Fourier Transform/Series.
Fourier analysis is of central importance in filtering and correlation because it allows efficient computation and gives theoretical insight into how these processes are related. The Fourier Transform is an analytical operation that allows a function that exists between the limits of – and + infinity to be expressed in terms of (complex) frequency and is a continuous function of frequency. However, a Fourier Transform of a real-World signal, which is sampled over a specific length of time – a record -, is not calculable. It can be approximated by a Fourier series, normally calculated through the discrete Fourier Transform algorithm, in which the signal is represented as the sum of a series of sine/cosine waves whose frequencies are an exact multiple of the fundamental frequency (=1./length of the record). Although this may seem to be splitting hairs, the differences between the Fourier Transform and the series are important. Fortunately, many of the relationships for continuous signals, for example a voltage wave form, are applicable to signals that are samples in time, which is the way that a signal is represented in a computer. The essential idea about the Fourier transform is that it takes a signal that is dependent on time, t, and represents it in a different domain, that of complex frequency, w. An operation performed on the time domain signal has an equivalent operation in the frequency domain. It is often simpler, far more efficient, and more informative to take a signal, convert it into the frequency domain, perform an operation in the frequency domain and then convert the result back into the time domain.
Some of these relationships are shown in figure 1.
Figure 1. The relationship between the input and output of a system in the time and frequency domains and their correlation functions. These are related through their (discrete) Fourier Transforms.
If a signal, x(t) passes through a linear system, typically a filter, the output, y(t) can be calculated from the input and the impulse response of the filter h(t), which is, mathematically, its response to an infinite amplitude spike that lasts for an infinitesimal time. The process by which this is calculated is called a convolution, which is often represented as “*”, so that:
y(t)=x(t)*h(t)
Looking at figure 1, this is shown in the blue upper panel. The symbol t, representing time, has a suffix, k, that indicates that this is a sampled signal at t0, t1, t2 …… Arrows represent the Discrete Fourier Transform (DFT) that convert the signal from the time to the frequency domain the inverse transform (DFT-1) that converts the signal from the frequency to the time domain. In the frequency domain, convolution is equivalent to multiplication. We can take a signal and transform it from x(tk) to X(wn). If we know the structure of the filter, we can calculate the DFT, H(wn), of its impulse response. We can write, using the symbol F as the forward transform and F-1 as the inverse transform, the following relationships to get the filter output:
X(w)=F[x(t)]
H(w)=F[h(t)]
Y(w)=X(w)H(w)
y(t)=F-1[Y(w)]
What we are doing is taking a specific frequency component of the input signal, modifying it by the frequency response of the filter to get the output at that frequency. The importance of the relationships shown above is that we can convert the frequency response of a filter, which is how filters are specified, into its effect on a period of a time domain signal, which is usually what we are interested in. These deceptively simple relationships allow the effects of a system on a signal to be calculated interchangeably in the time and frequency domains.
Looking at the lower panel in Figure 1, there is a relationship between the (discrete) Fourier Transform and the correlation functions of the inputs and outputs. The autocorrelation functions, which are the signals correlated with themselves are obtained by multiplying the transform by a modified form, the complex conjugate, written as X(w)*, (see below), which gives the signal power spectrum and taking the inverse transform. The cross correlation function is obtained by multiplying the DFT of the input by the complex conjugate of the output to obtain the cross-power spectrum, Gxy(w) and taking the inverse transform, i.e.:
Gxy(w)=X(w)Y*(w)= X(w)X*(w)H*(w)
Rxy(t)=F-1[Gxy(w)]
Therefore there is an intimate relationship between time domain signals representing the input and output of a system and the correlation functions of those signals. They are related through their (discrete) Fourier Transforms.
We now have to look in greater detail at the DFT, what we mean by a frequency component and what a cross-correlation function represents.
Figure 2 shows reconstruction of a waveform, shown in the bottom trace in bold by adding increasingly higher frequency components of its Fourier series. The black trace is the sum of the harmonics up to that point and the red trace is the cosine wave at each harmonic. It is clear that as the harmonics are summed, its value approaches the true waveform and when all the harmonics are used, the reconstruction is exact.
Up to this point, I have rather simplistically represented a Fourier component as being a cosine wave. If you compare harmonics 8 and 24 in figure 2, the peak of every third oscillation of harmonic 24 coincides with the peak in harmonic 8. In less contrived signals this does not generally occur.
Figure 2 A wave form, shown in bold, is constructed by summing its Fourier components shown in red. The black traces show the sum the number of Fourier components up to that harmonic.
The Importance of Phase
Each Fourier component has two values at each frequency. A sine and cosine waves are generated by a rotation of a point about an origin (Figure 3). If it starts on the y axis, the projection of the point is a sine wave and its projection on the x axis is a cosine wave. When the Fourier coefficients are calculated, the contribution of both a sine and a cosine wave to the signal at that frequency are determined. This gives two values, the amplitude of the cosine wave and its phase. The red point on the circle is at an arbitrary point and so its projection becomes a cosine wave that is shifted by a phase angle, usually written as j. Therefore the Fourier component at each frequency has two components, amplitude and phase and can be regarded as a vector.
Earlier, I glibly said that convolution is performed by multiplying the transform of the input by the transform of the impulse response (this is true since they are complex). This is equivalent to multiplying their amplitudes and adding their phases. In correlation, rather than multiplying X(w) and Y(w), we use Y(-w), the transform represented in negative frequency, the complex conjugate. This is equivalent to multiplying their amplitudes and subtracting their phases. Understanding the phase relationships between signals is essential in correlation[i].
Figure 3 The Fourier component is calculated as a sine and cosine coefficient, which may be converted to amplitude, A, and phase angle j. The DFT decomposes the time domain signal into amplitude and phases at each frequency component. The complex conjugate at is shown in blue.
Signal shape is critically determined by phase. Figure 4 shows two signals, an impulse and a square wave shown in black. I have taken their DFT, randomised the phases, while keeping their amplitudes the same, and then reconstructed the signals, shown in red.
Figure 4. The effect of phase manipulation. The black traces are the raw signal, and the red trace is the signal with a randomised phase spectrum but an identical amplitude spectrum
This demonstrates that phase has very important role in determining signal shape. There is a classical demonstration, which space doesn’t allow here, of taking broad band noise and imposing either the phase spectrum of a deterministic signal, while keeping the amplitude spectrum of the noise unaltered or doing the reverse: imposing the amplitude spectrum of the deterministic signal and keeping the phase of the noise unaltered. The modified spectrum is then inverse-transformed into the time domain. The phase manipulated signal has a high correlation with the deterministic signal, while the amplitude manipulated signal has a random correlation with deterministic signal, so underlining the importance of phase in determining signal shape.
Figure 5 The phase spectra of delayed impulses.
A very important concept is that phase represents delay in a signal. A pure time delay is a linear change in phase with frequency as shown in figure 5. The amplitude of the signal is unaltered, but in the case of a delay, there is increasing negative phase with frequency. However, any system that changes the shape of the input signal as it is passed through it, as is usually the case, will not have a linear phase spectrum. This is a particularly important concept when related to correlation.
We are now in a position to understand the cross-correlation function. Looking at the formula for correlation shown in figure 1, the CCF is:
,
which can be normalised to give the coefficient: 
, where N is the number of points in the record.
This rather formidable looking equation is actually quite straightforward. If k is zero, this is simply the standard formula for calculating the correlation coefficient and x is simply correlated with y. If k is one, the y signal is shifted by one sample and the process is repeated. We repeat this for a wide range of k. Therefore the function rxy is the correlation between signals two signals at different levels of shift, k and this tells one something about the relationship between the input and output.
We have a signal x(t) which has been passed through a physical system, with specific characteristics, which results in an output y(t) and we are trying to deduce the characteristics of the system, h(t). Since, from Figure1, the DFT, Y(w) of the output is the product of the DFTs of the input X(w) and the impulse response, H(w), could we not simply divide the DFT of the output by the DFT of the input to get the response? In principle, we can, providing the data is exact.
However most real world measurements contain noise, which is added to the inputs and outputs, or even worse other deterministic signals, and this renders the process somewhat error prone and the results of such a calculation are shown below (figure 6), illustrating the problem:
Figure 6. Left: Input (black) and output (red) signals for a system. Right: the calculated impulse response with 2.5% full scale amplitude noise added to the input and output (black) compared with the true impulse response (red). The low pass filtered response is shown in green.
This calculation illustrates another very important concept: linear physical systems store and dissipate energy. For example, a first order system, which could be the voltage output of a resistor/capacitor network or the displacement of a mechanical spring/damper system, absorbs energy transients and then releases the energy slowly, resulting in the negative exponential impulse response shown in figure 6. The variables which fully define the first order system are its gain and the time constant. The phase spectrum of the impulse response in distinctly non-linear. Attempts to measure another variable, for example delay, which implies a linear phase response, is doomed to failure because it doesn’t really mean anything. For example, if one is looking at the relationship between CO2 and temperature, this is likely to be a complex process that is not defined by delay alone and therefore the response of the system should be identified rather than a physically meaningless variable.
Noise and Correlation
Correlation techniques are used to reduce the effects of noise in the signal. They depend on the noise being independent of, and uncorrelated with, the underlying signal. As explained above, correlation is performed by shifting the signal in time, multiplying the signal by itself (auto-correlation) or with another signal (cross-correlation), summing the result, and performing this at every possible value of shift.
Figure 7. Broad band noise ( black) with the autocorrelation function( red) superimposed.
In broad band noise, each point is, by definition uncorrelated with its neighbours. Therefore, in the auto-correlation function, when there is no shift, there will be perfect correlation between it and its non-shifted self. For all other values of shift, the correlation is, ideally, zero, as shown in figure 7.
The auto correlation function of a cosine wave is obtained in the same manner. When it is unshifted, there will be a perfect match and the correlation will be 1. When shifted by ¼ of its period, the correlation will be zero, be -1 when shifted by ½ a period and zero when shifted by ¾ of period.
The ACF of a cosine wave is a cosine wave of the same frequency with an amplitude of the square of the original wave. However if there is noise in the signal, the value of the correlation will be reduced.
Figure 8 shows the ACF of broadband noise with two sine wave embedded in it. This indicates recovery of two deterministic signals that have serial correlation and are not correlated with the noise. This is a basis for spectral identification in the presence of noise.
Figure 8 The ACF of a signal containing two sine waves of different frequencies embedded in noise. The ACF (red) is the sum of two cosine waves with the same frequencies.
A very important feature of the ACF is that if destroys phase information. Referring to Figure 1, the DFT of the ACF is X(w) (or Y(w)) multiplied by its complex conjugate, which has the same amplitude and negative phase. Thus when they are multiplied together, the amplitudes are squared and the phases are added together making the resultant phase zero. This is the “power spectrum” and is the ACF is its inverse DFT. Therefore the ACF is composed entirely of cosine waves and is symmetrical about a shift of zero.
However, the cross-correlation function, which is the inverse DFT of the cross-power spectrum contains phase. By multiplying the complex conjugate of the output by the input in the frequency domain, one is extracting the phase difference and the delays at each frequency between the input and the output and the cross-correlation function reflects this relationship. If the power spectrum, e.g.: the DFT, of rxx(t) is Gxx(w) and the cross-power spectrum of rxy(t) is Gxy(w), then:
H(w)=Gxy(w)/Gxx(w)
Figure 9 shows the same data as used in figure 6 to calculate the impulse response and the error is very much reduced because signal correlation is a procedure that separates the signal from noise.
Figure 9. Estimate of the impulse response using the data in figure 6 via cross-correlation (black) and pre-filtering the data (green).
These calculations are based on a single record of the input and output. When available, one uses multiple records and calculates the estimated response E[H(w)] from the averaged power spectra:
E[H(w)]=< Gxy(w)>/<Gxx(w)>
Where < x> means the average of x. This leads to a better estimate of the impulse response. It is possible to average because correlation changes the variable from the time domain to relative shift between the signals so aligning them.
One simple check that can be performed to check that one is getting reasonable results, assuming that one has enough individual records, is to calculate the coherence spectrum. This is effectively the correlation between the input and output at each frequency component in the spectrum. If this is significantly less than 1.0, it is likely that there is another input, which hasn’t been represented, or the system is non-linear.
One of the major problems in applying signal processing methods to climate data is that there is only one, relatively short, record and therefore averaging cannot be applied to improve estimates.
Improving resolution by record segmentation and filtering.
One can improve estimates of the response if one has a model of what the signal represents. If one is dealing with a short term process, in other words the output varies quickly and transiently in response to a short term change in input, one can estimate the length of record that is required to capture that response and hence the frequency range of interest. This enables one to segment the record into shorter sections. Each segment has the same sampling frequency, therefore the highest frequency is preserved. By shortening the length of each record we have thrown away low frequencies because the lowest frequency is 1/(record length). However, we have created more records containing high frequencies, which can be averaged to obtain a better estimate of the response.
The other strategy is filtering. This, again, involves assumptions about the nature of signal. Figure 10 shows the same data as in figures 7 & 8 after low-pass filtering. The ACF is no longer an impulse but is expanded about t=0. However the ACF of the deterministic signal is recovered with higher accuracy.
This can be done here because the signal in question has a very small, low frequency, bandwidth and is not affected by the filtering (figure 11). The effects of the filter are easily calculable. If it has a frequency response of A(w), the input and output signals become X(w)A(w) and Y(w)A(w). The cross correlation spectrum is therefore simply:
Gxy(w)= A(w)X(w)A*(w)Y*(w)=A2(w)X(w)Y*(w)
Figure 10 The ACFs using the same data as in Figure 6. Note that the ACF of noise is no longer an impulse at t=0 and that there has been a considerable improvement in the ACF of the two sine waves as they now represent a higher fraction of the total power in the signal.
A2(w)is the autocorrelationfunction of the filter, which has no phase shift and will not affect the phase of the cross-power spectrum, provided the same filter is used on the input and output. This is because the phase reversal of the complex conjugate of the filter in the output cancels out that applied to the input, so the timing relationships between the input and output will not be affected. Provided the ACF spectrum is in the pass band of the filter, it will be preserved. In figure 9, the estimated impulse responses are shown using filtered (green) and non-filtered data. If one wishes to characterise the response, by assuming it is a first order system (which this is), one can fit an exponential to the data so getting its gain and time constant. The filtered result gives marginally better estimates but one has to design the filter rather carefully, appreciate that filtering modifies the impulse response and correct the results for this.
Thus, it is possible to filter signals and perform correlations, provided that the frequency band being filtered does not overlap to much the system response, as illustrated in figure 11, and one is careful.
Figure 11. The signal spectrum is composed of the true signal and noise spectra. A good filter response (green) will attenuate (grey) some of the noise but preserve the true signal, while a bad filter will modify the true signal spectrum and hence the ACF.
In practice, however there is likely to be an overlap between the noise and signal spectra. If the system response is filtered, the correlation functions will be filtered and become widened and oscillatory. In this case, the results won’t mean much and almost certainly will not be what you think they are! There are more advanced statistical methods of determining which part of the spectra contain deterministic signal but, in the case of climate data, the short length of the modern record and the poor quality of the historical record makes this very difficult.
Degrees of Freedom.
Suppose we have two signals and we want to determine if they have different means. They both have a normal distribution and the same variance. Can we test the difference in means by performing a “Student’s t” test? This will almost certainly be wrong, because in most simple statistical tests, there is an assumption that each observation is independent. In figure 7, the ACF is an impulse and nominally zero elsewhere, showing that each point is independent of each other. If the signal is filtered, the points are no longer independent because we have convolved the signal with the impulse response of the filter, as shown in figure 10. Looking at figure 1, the time domain convolution is given by:
This is similar to the correlation formula, except that the impulse response is reversed in time. It shows that the output at any point is a weighted sum of the inputs that have preceded it and are therefore no longer independent. Therefore in applying statistical tests to signal data, one has to measure the dependence of each sample on others by using the autocorrelation of the signal to calculate the number of independent samples or “degrees of freedom”.
Conclusion.
Correlation and filtering are highly interdependent through a set of mathematical relationships. The application of these principles is often limited because of signal quality and the “art” of signal processing is to try to get the best understanding of a physical system in the light of these constraints. The examples shown here are very simple, giving well defined results but real world signal processing may be messier, require much more statistical characterisation and give results that may be limited statistically by inadequate data.
One always has to ask what is the goal of processing a signal and does this make any sense physically? For example, as discussed earlier, cross correlation is often used interchangeably with “delay and it is only meaningful if one believes that phase response of the system in question has a linear phase response with frequency. If one is estimating something that is not meaningful, additional signal processing will not be helpful.
Rather, if one has a model of the system, one can then calculate the parameters of the model and, having done this, one should look carefully at the model to see if it accounts for the measured signals. Ideally, this should be tested with fresh data if it is available, or one can segment the raw data and use one half to create the model and test it with the remaining data.
Modern scripting programs such as “R” allow one to perform many signal processing calculations very easily. The use of these programs lies in not applying them blindly to data but in deciding how to use them appropriately. Speaking from bitter experience, it is very easy to make mistakes in signal processing and it is difficult to recognise them. These mistakes fall into three categories, programming errors, procedural errors in handling the signals and not understanding the theory as well as one should. While modern scripting languages are robust, and may largely eliminate straight programming errors, they most certainly do not protect one from making the others!
[i] I have used the electrical engineering convention in this post, i.e.: the basis function of the Fourier Transform is a negative complex exponential.
Related articles
- Fourier series before Fourier (johndcook.com)
- Retrotechtacular: The Fourier Series (hackaday.com)
- Sampling Rate (tamino.wordpress.com)
RCSaumarez says:
April 10, 2013 at 2:45 pm
” In which case you will end up with a practically random signal. “
How do you get there? It isn’t random at all. The SNR is high. Look at the plot.
Do you even know what SNR is?
I GOT my PhD in this stuff, and you need schooling.
Greg Goodman says:
April 10, 2013 at 2:50 pm
“…so I guess that’s a no.”
You mean this plot? And, we are talking about this relationship? Wow, that Bart guy and I have a lot in common.
RCSaumarez says:
April 10, 2013 at 2:51 pm
God, how I had smug ignorance. Dude, you have no clue.
RCSaumarez says:
April 10, 2013 at 2:51 pm
God, how I hate smug ignorance.
No, it was not Freudian. This is outrageous.
“You appear to be unable to do elementary calculus and your derivation of a the frequency response of a differebtiator is just plain wrong”
What the hell do you THINK the derivative operator is for a Laplace transform???
Here’s a clue: d/dt(exp(s*t) = s*exp(s*t).
What kind of an idiot makes a statement like yours with no understanding whatsoever???
Why don’t you at least look it up in a few links before spouting off so stupidly?
“You mean this plot? And, we are talking about this relationship? Wow, that Bart guy and I have a lot in common.”
you flatter yourself.
RCSaumarez says:
April 10, 2013 at 2:51 pm
Maybe you are confused about the discrete time differentiation operator, D(z) = (1 – z^-1)/T. This is what the WTF site is doing when it “differentiates”.
I doubt you have a clue what the Z-transform is. I will walk you through it.
With z = exp(s*T), get
D(exp(s*T)) = exp(-s*T/2)*(exp(s*T/2) – exp(-s*T/2))/T
Evaluated at s = j*w, get
D(exp(j*w*T)) = 2*j*exp(-j*w*T/2)*sin(w*T/2)/T
In order to line up the data in time, you have to advance the phase (shift the time series forward 1/2 step). Thus, the response for the phase advanced discrete time derivative (which is what the WTF site implements) is
D(exp(j*w*T))*exp(j*w*T/2) = 2*j*sin(w*T/2)/T
THERE. Do you understand NOW???
Greg Goodman says:
April 10, 2013 at 2:50 pm
“OK, so eye-balling the graph is all you’ve got.”
You say this as if it were a bad thing. Quite the contrary, when you can look at data and see an obvious relationship, that is a very good thing.
Whatever method of extracting some numerical figure of merit you might come up with, it will necessarily reflect that obvious relationship. Indeed, if it did not, you would know from that sanity check that you had done your analysis wrong.
Everything I have stated is so elementary and obvious, it just leaves me nonplussed to step into this bizarro, surreal world and be assualted with arguments such as I have sustained on this thread. In my world, among my peers, I am accounted an authority. In this forum, I am the proverbial one eyed man in the land of the blind who, contrary to the usual saying, is no king, but a lunatic spouting off that he can “see” things, whatever that means.
Well, anyway, we’ve gone about as far as we can go. Good luck in your future endeavors, Mssrs. Goodman and Saumarez.
You’re going to need it.
Greg Goodman says:
April 10, 2013 at 2:50 pm
“OK, so eye-balling the graph is all you’ve got.”
One last thing. As I have pointed out on numerous occasions above:
dCO2/dt = k*(T – To)
The gas is continuously upwelling from ocean circulation, volcanism, and biological process. It is continually downwelling from ocean circulation, mineral weathering, and biological processes. And, it is doing so at some rate. Any differential in rate between what is coming in, and what is going out, is either going to result in either an increase or decrease in the surface system.
The rates at which it comes in, and goes out, are strongly affected by temperatures, and not generally in the same measure.
So, you have a temperature modulate rate of increase or decrease. That is essentially what the equation says: the rate of change is proportional to the temperature anomaly relative to a particular baseline where equilibrium attains, i.e.,
dCO2/dt = k*(T – To)
And, lo and behold, we look at the data, and that is exactly what we see.
Come on, guys. This is unbelievably straightforward. How can you not see it?
Bart is absolutely correct. Seasonal temperatures activate CO2 and anthropogenic inputs of CO2 fill in the overall secular trend.
No Bart, I don’t understand, or rather I do.
f=sin(wt), let u=wt,
df/dt =d(sin(u))/dt . du/dt= w .cos(wt).
This is an increase in amplitude of w and a phase shift of pi/2.
AS regards my future endeavors, there is quite a large body of opinion that I have been reasonably successful over the last 20 years.
sorry, that should read:
df/dt=d(sun(u))/du. du/dt = w cos(wt)
RCSaumarez says:
April 11, 2013 at 4:18 am
“AS regards my future endeavors, there is quite a large body of opinion that I have been reasonably successful over the last 20 years.”
As have I. In, I might add, THIS particular field. Prudence, if not courtesy, should dictate that you take that into consideration before jumping to conclusions.
You seem to have a confusion regarding continuous time differentiation versus discrete time first differencing and dividing by sample time.
Do you think the WTF site is performing continuous time differentiation on its data? How, precisely, would it accomplish this? How are you going to get amplification beyond the Nyquist frequency for discrete data?
Auto Phil says:
April 11, 2013 at 3:59 am
“Seasonal temperatures activate CO2 and anthropogenic inputs of CO2 fill in the overall secular trend.”
The trend is fully accounted for by the temperature. Anthropogenic inputs are superfluous. They are unnecessary for the reconstruction of CO2 levels.
Look, guys, this match is signal, not noise. It is not random. These are two entirely different time series, and yet they lie right on top of each other. The odds of that happening randomly are effectively zero.
Both the trend, and the variation, match as perfectly as likely possible with inherent measurement errors. There is no room to add in additional forcing from anthropogenic inputs without subtracting out the long term temperature trend. And, if you do that, the variation will no longer match. Occam’s Razor: temperature is driving essentially the whole shooting match, and anthropogenic inputs have, at most, a minor role.
RCSaumarez says:
April 11, 2013 at 4:18 am
f=sin(wt), let u=wt,
df/dt =d(sin(u))/dt . df/dt= w .cos(wt).
Let f = sin(w*T). The phase advanced discrete time differentiation is
df/dt = ( f(t+T/2) – f(t-T/2) ) / T
= (sin(w*(t+T/2)) – sin(w*(t-T/2)) / T
= cos(w*t) * ( 2 * sin(w*T/2) / T)
Let f = sin(w*t) in the above.
RCSaumarez says:
April 10, 2013 at 2:45 pm
“If I understand you correctly. which some might find difficult, you have differentiated a low pass fileered signal. In which case you have attenuated the deternistic components and amplified the noise. In which case you will end up with a practically random signal.”
Looking back up the thread, I realized from this statement that you totally misapprehended what is going on. By low pass filtering, one has not “attenuated the deternistic components and amplified the noise.” One has attenuated the noise and passed through the deterministic components. Deterministic components are generally low frequency, and noise is generally high frequency.
Bart,
I’m beginning to see where you are coming from.
The problem with phase advanced differentiation, is that it is not a true derivative. It is a band pass filter with zeros at the fundamental and at the Nyquist frequencies and a pi/2 phase shift across the frequency range. In the lower frequency ranges, it is an approximate derivative in terms of amplitude, except that it does not have the phase characteristics of a derivative.
That’s OK unless you want to treat the results as the driving function of a differential equation. In this case you a substituting some in place of a derivative, when it is definitely not.
In your differential equation:
dCO2/dt=k(t-t0)
what is the significance of k? Written like this I would have assumed that k is a constant. Do you intend it to be an arbitrary (complex) function, in which case, this equation has a very different complexion.
If k is a constant, why not simply integrate temperature and correlate it with CO2?
Looking at the plot of dCO2/dt and temperature in the link, there is certainly a relationship between them. However, it is also clear that there is a variable lead/lag between them. Given the solubity of CO2, might one not expect ther to be a non-linear relation between them? In which case, any simple linear processing is not useful.
On a wider point, I really can’t see the point you are trying to make.
I’m sorry, I made a mistake typing, the center difference derivative DOES have the phase characteristics of a differentiatior, i.e. a pi/2 phase shift but over the frequency band only approximates in the lower frequency range. This was a slip of the keyboard. Apologies.
We should applaud Bart for finding a new chemistry whereby a small increase in temperature will effectively drive all the CO2 out of the ocean. That is such a neat trick.
Bart. If you really want to characterise the relationship, and since the CO2 is a trend, have you considered a state variable approach, i.e.: Kalman Filters?
Friends:
Richard Saumarez provided a fine article which gives a good overview of Fourier analysaes. So, I have been following your discussion with interest in genuine hope of learning more about principles of time series analyses.
I write to respectfully suggest that you are going in a tangent which is not helpful to those – including me – observing this thread with a view to learning. And I am writing in hope the discussion will get back ‘on track’.
Before explaining my point, I need to make something clear. Bart knows I do not accept the indications of his analysis, but that is NOT why I am asking for the thread to stop being side-tracked onto debate of his analysis. Indeed, I do not like snark about Bart inventing a new chemistry: if it is right then his analysis indicates biological changes probably in the ocean surface layer.
Detailed consideration of a case would be a useful expansion of the article by Richard Saumarez, but Bart’s analysis is far too limited a case to provide such expansion.
If you want to use investigations of temperature and CO2 time series as a case study then it would be useful to onlookers – such as me – to compare various ways of conducting such analyses, and not the merits of only one.
However, temperature and CO2 may not be the optimum such case study for illustration of the principles outlined in the article by Richard Saumarez. Only Richard Saumarez can decide that.
I hope everybody will understand I am writing this with sincerity and my only intent is to assist the value – now and as an archive – of this thread.
Richard
@ur momisugly Richard S Courtney
Thanks, I felt that the whole drift of the conversation was going off the rails. What started with a rather simple minded question on whether it kosher to filter a signal and then form a correlation has expanded well beyond its remit into something that is far more tricky.
I am not sure that signal processing for its own sake is useful here. There are (at least ) two issues.
The first is whether temperature drives CO2 (or its derivative) or vice versa. Is the system linear, or at least linearisable? Does the basic physical chemistry and physics suggest that it is? Should, as you suggest consider that there is a biological compartment that absorbs CO2 rapidly and releases it much more slowly, if at all? Therefore, without a model, I am not sure that DSP would be helpful. I’m no expert in this field – I do biomedical engineering.
The second issue is the length and quality of data. Although modern data is better than the historical record, one is dealing with a short segment segment with a lot a variability. Normally, one forms ensembles (i.e. averages in English) to improve estimates. Looking at the Temperature CO2 data, there is clearly some relationship, but there is also considerable leads and lags between them and any estimate would have large error limits in determining this – I’m guessing but I suspect that one wouldn’t get an answer that was helpful. In my post I remarked that one of the major problems in climate based signal processing is that we only have one, relatively short record which precludes the normal approach where one tries to gets tons of data, divide it into many seperate records whose length captures the effect one is looking for, so that one can do sensible statistics. ( Rather counter-intuitively having a long record does not increase the accuracy of estimates. It simply enables one to get an single estimate over a wider frequency range. The key is being able to average data). There are undoubtly technical issues on how you handle the signals, but these are details – although it is important that these are done correctly.
I would say that signal processing is a tool that allows one to test hypotheses about the structure of a system rather than a rod and line that allows one to go on a fishing expedition. Once one has a hypothesis that is based on sensible assumptions, one can make statements about the quality of the data required to test that hypothesis. My feeling is that the data isn’t good enough and that any simple model won’t account for their relationship..
Auto Phil says:
April 11, 2013 at 2:57 pm
“We should applaud Bart for finding a new chemistry whereby a small increase in temperature will effectively drive all the CO2 out of the ocean.”
There is A LOT of CO2 in the ocean. A small differential rate between upwelling and downwelling, integrated over time, makes a huge difference to the surface system.
RCSaumarez says:
April 11, 2013 at 2:31 pm
“…what is the significance of k?”
It may be taken to be constant over a given interval of time. It holds up pretty well as a constant since 1958, when precise CO2 measurements began to be made. But, overall, it should probably vary over time, possibly not even particularly smoothly.
For example, in the circulation of the oceans, there may be different CO2 concentrations along the flow in the “pipeline”. When a relatively (to overall surface waters) CO2 rich volume starts surfacing, it would start outgassing and accumulating in the atmosphere. The constant “k” would increase, and the effective equilibrium temperature To would shift downward, as temperatures would have to cool significantly to allow equilibrium between inflow and outflow to be reestablished.
But, all these dynamics would be slow, and within the space of some years, the parameters could be considered relatively constant. The simple model
dCO2/dt = k*(T – To)
is then seen as a linearization of a more complicated, nonlinear and time varying process.
The units of “k” would be ppm/degC/unit-of-time, and basically is simply the amount by which the rate of accumulation in the atmosphere changes per degree of temperature change.
“If k is a constant, why not simply integrate temperature and correlate it with CO2?”
You could do that, but as an integrator attenuates high frequencies, you will miss a lot of detail, and that is the discriminator which tells us that this is really what is driving atmospheric CO2, and not human inputs.
If you look at anthropogenic rate of input, you will see it does not really correlate with the rate of change of atmospheric CO2. The total accumulation of anthropogenic inputs is slightly quadratic, and so also looks affinely similar to the measured CO2 concentration. But, it is only a superficial resemblance, as the rate of change does not correlate to the measured rate of change like the temperature does.
The superficial resemblance of accumulated anthropogenic inputs to measured concentration is really not remarkable or compelling. Both are slightly quadratic series, and you can always make two slightly quadratic series look affinely similar simply by doing a linear regression of the one against the other to find the affine parameters which produce the best fit. It is when you get into the fine detail of the rate of change that you recognize that the resemblance is superficial.
“However, it is also clear that there is a variable lead/lag between them. “
Keep in mind that the WFT site is doing the numerical differentiation and averaging, and advancing the resulting time series to match up in phase. So, it is an anti-causal filter. Therefore, events in the CO2 series can appear to have started before the time at which they actually do.
Moreover, it is important to keep in mind that the temperature measure is a global spatial average, and whatever events are actually causing the CO2 rise may not be in evidence in that average until after the CO2 measurement has been taken. Some people have noted that there appears to be an especially good fit with Southern Hemispheric data, which suggests this may be largely an oceanic phenomenon.
“On a wider point, I really can’t see the point you are trying to make.”
If you do your correlation analysis between the rate of change of CO2 and the temperature, you will find that the dynamics are a good fit to the dCO2/dt = k*(T – To) model, and you may discover finer structure which could lead to better understanding of how this comes about. I don’t do this stuff for a living, and do not have time to chase down everything. Moreover, I am sure I am not nearly as familiar with the particular system as you are, and you can most likely come up with connections of which I am not aware which would help explain how this all comes about.
RCSaumarez says:
April 11, 2013 at 2:59 pm
“If you really want to characterise the relationship, and since the CO2 is a trend, have you considered a state variable approach, i.e.: Kalman Filters?”
Absolutely. But, I do not have the time. I wish someone would do it.
richardscourtney says:
April 11, 2013 at 3:30 pm
My goal was not to sidetrack the conversation. Quite the contrary, I believe the CO2/temperature relationship is an excellent case study for the methods Richard has outlined, and that is why I brought it up.
RCSaumarez says:
April 11, 2013 at 4:58 pm
“The first is whether temperature drives CO2 (or its derivative) or vice versa.”
If the rate of CO2 drives temperature, then it can accumulate to very high concentration but, once it stops accumulating, the temperature will drop back to the equilibrium level.
That does not make any sense. Thus, we are compelled to conclude the converse.
richardscourtney says:
April 11, 2013 at 3:30 pm
“Bart knows I do not accept the indications of his analysis…”
I still have no idea how you can imagine that this kind of correlation between two completely independently measured quantities can arise by chance, and not be a real, shared characteristic.
Bart says:
April 11, 2013 at 5:04 pm
“…once it stops accumulating, the temperature will drop back to the equilibrium level.”
Even more absurd, you could ramp the CO2 up an insignificant amount, but in a very short time, and temperatures would spike upward.
No, the rate of CO2 accumulation is not driving temperature. Temperature is driving the rate of change of CO2. Allow me to remind everyone of the obvious dynamic. From a previous post:
BTW, I have what I believe to be solid reasons to suspect that a fuller description of this system is something along the lines of
dCO2/dt = (Co2eq – CO2)/tau + f*H
dCO2eq/dt = k*(T – To)
where CO2eq is the equilibrium level of CO2 to which the system is attracted based on current conditions, tau is a time constant associated with sequestration, H is the rate of human inputs, and f is the fraction of human inputs which is not rapidly absorbed into the oceans and remains airborne. If tau is “short”, then H will be rapidly sequestered, and CO2 will track CO2eq.
The quantities “f” and “tau” could be made operator theoretic to account for more complex dynamics, but I think starting with an assumption of constant values and attempting to estimate them from the data would be a good starting point.
Bart is right.
His term H is the human forcing function for CO2 and it swamps the temperature activated portion.