A faulty tapering

I see there's renewed interest at Climate Audit in discussing Bart's code for relating temperature and cloud radiative fluxes (CRF). I've discussed it a lot there, and on the past few posts here. But I'm bogged down there explaining what I believe is a simple error, which I'd like to explain here in more detail.

There is questioning at CA as to whether the whole approach is the right way to try to analyse feedback between these variables. There is a confused notion of causality. I agree with that. But here the aim is more limited.



Here's the problem. We have a variable dR (CRF) and we have a model for its dependence on T: \[dR(t) = \int_{-\infty}^\infty h(t-\tau) T(\tau) d\tau\] Just convolution. This can be resolved by Fourier transform FT, which converts convolution into a product: \[FT(dR) = FT(h) * FT(T)\] and so \[h = iFT(FT(dR)/FT(T))\] iFT being the inverse FT. Of course, all this is implemented wqith the discrete Fourier transform (DFT), and in practice with FFT.

I'll analyse a simplified version of Bart's code translated to R. Here it is:

data=read.csv("flux.csv",skip=0);

data[is.na(data)]=0;

temp=data[,9];

dR=data[,5]-data[,8];

# temp and dR are the data vectors, length 124 months N=length(dR);

dT=1/12;

Nsamp = 8192;

Npad=Nsamp-N;

X=fft(c(temp,rep(0,Npad))); # pad temp to 8192 points and FFt

Y=fft(c(dR,rep(0,Npad)));

# Inverse FFT to get the impulse response
h = Re(fft(Y/X,inv=T))/dT/Nsamp;

# now h is tapered; Nc=Nsamp/8;

w=c(rep(1,Nc),(1+cos(pi*(1:Nc-1)/(Nc-1)))/2,rep(0,6*Nc));

H=dT*fft(h*w); # Fourier transform of tapered h;

So X and Y are duly FFT'd, with lots of padding. There are 8192 frequencies, all multiples of a base frequency (1/8192) RPM (revs per month), and 8192 month in the time range.

Periodicity

Now the first thing to remember with the DFT is that everything is periodic, with period 8192 units. That is because everything is expressed in linear combinations of harmonics of that base frequency. That applies to everything that emerges from an FFT or an iFFT. Being periodic, it is best to visualize it on a circle.

Impulse response

So back to the code. Down to the calculation of h, as an FFT of Y/X, it's straightforward. And h is indeed the vector which convolves with T to give dR to near machine accuracy. To emphasise the periodicity, I'll show two periods. It looks like this:


The y-axis is scaled for the smoothed h (30 month triangular filter). The unsmoothed is scaled down by a factor of 40. Now I'll show a window about t=0:


There's clearly a lot of noise. The reason is basically that there are only 104 units of real information of temp, and the large number of vectors orthogonal to that are used to fit the padding.

Both sides matter

Now Bart says that the negative t part of h is pure noise and can be discarded. That's not true. To show this, just run the problem with the data in reverse order. The problem has exactly the same form. What does the smoothed h then look like:


The same as before, but exactly reversed. Is the negative t part now noise? It's exactly the same numbers as before, on the positive side.

Very strange taper

So Bart wants to taper h. He creates a taper w which looks like this. Again it's periodic (there's no other domain), so I'll show 2 periods:


Max w is 1; the y axis is scaled for h, which is shown in the background.

So now there are two problems. One is that the taper cuts of a part of h that matters - the negative t part. This has the result that the tapered h no longer convolves with T to give dR. In fact, it gives this, after a little smoothing (10 mth triangular):



The red curve from h exactly matches; the gold curve from tapered h does not. That's serious, because reproducing dR is what h is supposed to do. If it reproduces something else, then any results apply to that quantity.

Frequency response of taper

A less serious, but even more obvious issue is the taper itself. When you cut data, you in effect multiply by a gate function, which is 0 where you cut, and 1 elsewhere. The effect in the frequency domain is to convolve with the FT of the gate. You want to make as little difference as possible in medium frequency ranges, so that FT should be narrow, like a delta function. But an actual gate (square wave) tapers off slowly in the frequency domain, as inverse freq, often described as 6 dB/octave.

The rule is that the rate of roll-off (in dB) is proportional to the order of the first derivative that does not exist (somewhere). For a gate, that is 1. For a Hann taper, which is a  period of cos^2, the second derivative exists, but is discontinuous at the join. So the roll-off is 18 dB/octave.

Now Bart is using half a Hann filter, with an adjoining flat stretch. That's fine, but is spoilt by the cut at zero. In the frequency domain, his taper looks like this (magnitude). The red line has slope 6 dB/octave


As a taper, it performs very poorly. It should attenuate at at least 18dB/octave. Now of course that is not so serious. It could be omitted. But it indicates the wrongness of the thinking.

For completeness, here is the R graphics code that I used:

# For display, I stretch the variables to cover two periods.

# Then I select subsections of the range using xlim

f1=c(1:15,15:1); f1=f1/sum(f1);

q=Mod(fft(w)); q=c(q,q);

w=c(w,w);h=c(h,h);

hs=filter(h,f1,circ=T);

hm=mean(h[1024:7168]); hh=h-hm;

f1=c(1:5,5:1); f1=f1/sum(f1);

cum=cumsum(hh)*dT; cum=filter(cum-cum[Nsamp/2],f1);

x=(1-Nsamp):Nsamp;

wh=w*h; wm=w*(max(wh)/2)

xl=c(1,Nsamp,-300,300,-800,800);

graphics.off();

cl=c("black","grey","red","green")

png("post/f1.png",width=1000);

plot(x,h/40,type="l",col=cl[2],ylim=c(-2.5,1.5),xlab="Month",ylab="Response W/m2/C",

main="Impulse response CFR to temp");

lines(x,hs,col=cl[3]);

legend("topright",c("h","h smoothed"),text.col=cl[2:3])

dev.off()

png("post/f2.png",width=1000);

plot(x,h/40,type="l",xlim=xl[3:4],ylim=c(-2.5,1.5),col=cl[2],xlab="Month",ylab="Response W/m2/C",

main="Impulse response CFR to temp");

lines(x,hs,col=cl[3]);

legend("topright",c("h","h smoothed"),text.col=cl[2:3])

dev.off()

png("post/f3.png", width=1000);

plot(x,wh,type="l",col=cl[2],xlab="Month",ylab="Response W/m2/C",

main="Tapered impulse response CFR to temp, and taper");

lines(x,wm,col=cl[3]);

legend("topright",c("Tapered h","taper"),text.col=cl[2:3])

dev.off()

png("post/f5.png");
plot(x/Nsamp,q,type="l",xlim=c(1,800)/Nsamp,log="xy",col=cl[1],xlab="Freq RPM",ylab="Mag FFT(taper)",
main="Impulse response CFR to temp");
lines(x/Nsamp,1000/x,col=cl[3]);
dev.off()