Monday, September 12, 2011

Impulse responses and the Spencer/Dessler Cloud models

I've been involved in discussions at this CA thread, which has gone into FFT analysis of the relation between temperature and cloud contribution to TOA energy flux, as discussed a lot re papers from Dessler and Spencer and Braswell.

I should hasten to add that the analysis at CA is based on the use of all CERES data, which Dessler did not do. He gave good reasons for his choice to use reanalysis instead. Consequently, I don't think this alternative approach is telling us anything useful about clouds. However, there is interesting maths.

Commentator Bart there proposed a model in which an impulse function is found by FFT which, when convolved with surface temp (Hadcrut3) will reproduce that "cloud forcing" function, in his notation dR (for δR_cloud). He has produced a lot of Matlab analysis. His contention is that the low frequency behaviour of the impulse response shows strong negative feedback. I don't now agree with that - I think he is looking at very low frequency results which can't be supported by the short period (about 10 years) of data. Furthermore the integral of the impulse response that he cites is the low frequency limit where it becomes just the ratio of the time gradients of the two data sets as determined by OLS regression.

Anyway, I translated his code into R and did some comparison calculations, which I put n a page attached to this blog. I'm now transferring it to this post, in case people want to comment. I'm hoping to write a more detailed analysis soon. The data being used, referred to in the code as flux.csv, is here.



I have converted Bart's matlab code to R - see below. I've only done as far as the impulse response so far.
I'm testing a Hanning window - a cos^2 function which tapers to zero at both ends of the data.

When you FFT a finite set of data, it's as if you did the full set with a gate function. You have to accept low frequency problems with finite data, but the gate adds HF signal as well. A Hanning window tapers much faster in the freq domain, and mitigates this extra problem.

First a check that my R code is getting Bart's result. I've commented out the taper line:





Now after applying the Hanning taper to temp and dR

Update. Steve at CA asked to see if these impulse responses could be convolved with the temp to give dR, which is really what they are for. It's a useful check. In R, you have to use type="circular" and conj=F in the convolve() options. That done, it works. In the plots below, black is the original dR (cloud radiance), red uses the unsmoothed impulse response, and the cyan is the smoothed impulse response. I made the black wider so you can see the red superimposed.

The plot on the right is the same, but smoothed with a 10 month triangle filter. This was done without the Hann tapering of the data. I should clarify that the smoothed impulse response referred to here is the top graph. I have not showed the unsmoothed response - perhaps Bart has. 

Update I've added in dark gold the result of convolving with the taper (w) that Bart used on h. This taper makes h effectively one-sided in time.




Here's the R code. I've tried to follow Bart's notation and steps.
Update - the code originally here deviated in minor ways - I've added fixes

data=read.csv("flux.csv",skip=0);
data[is.na(data)]=0;
temp=data[,9];
dR=data[,5]-data[,8];
N=length(dR);
h=abs(1:N-64.5)/62.5
h[h>1]=1;Ha=cos(h*(pi/2))^2;  #  Hanning Window
temp=temp*Ha; dR=dR*Ha; #  Tapering
dT=1/12;
Nsamp=8192/2;
Npad=Nsamp-N;
X=fft(c(temp,rep(0,Npad)))+1.0e-9;
Y=fft(c(dR,rep(0,Npad)));
h = Re(fft(Y/X,inv=T))/dT/Nsamp;
Nc=Nsamp/8;
w=c(rep(1,Nc),(1-cos(pi*(1:Nc-1)/(Nc-1)/2)),rep(0,Nsamp-6*Nc))
hw=h*w;
f1=c(1:15,15:1); f1=f1/sum(f1);
hs=filter(h,f1)
t=(0:599)*dT
png("impresHanning.png");
plot(t,hs[1:600],type="l",ylab="Impulse Response W/m2/C/yr",xlab="Years",main="Cloud-Temperature System Smoothed Impulse Response/Hanning")
dev.off()


Late final extra

Carrick asked to see the mag plot of h. I've followed Bart's scheme, with log axes. It's weighted towards low frequencies. This h has neither my Hann taper of data nor Bart's h truncation. Mag and phase:



10 comments:

  1. Could you show what the amplitude of "h" (the transfer function) looks like?

    Also, it's pretty conventional when doing this sort of thing to look at the coherence between the two series. It'd be interesting to see that displayed too, if you have time.

    ReplyDelete
  2. Carrick,
    I've put them at the bottom of the post.

    I agree cross-spectral is probably more informative. I'll try to get that working.

    ReplyDelete
  3. Hi Nick,
    I'm confused here. When you put up the seperate html page initially, you got a replication of Barts smoothed impulse response which looked very similar to Bart's result
    http://climateaudit.files.wordpress.com/2011/09/cloudimp.jpg

    But now you seem to have removed the smoothing so it looks a mess compared to your version with the Smoothed/Hanning filter.

    Why have you done that?

    ReplyDelete
  4. TB,
    No, it hasn't really changed, though I should have explained it better. What I've done is to include the negative time material. This is a bit controversial. But if you look at the zero on the X axis and what is to the right, it's very similar to what Bart shows. The small difference is right near x=0; the reason for that is that I've smoothed across zero. Bart's curve fades out at that point, because he doesn't regard the values at zero and below as meaningful.

    You can see the similarity in the next diagrams down, where that LH stuff is removed.

    However I admit there is a problem. I changed the plits for the next post, and explained it there, but didn't note it here.

    I've stopped (for the moment) using the Hann filter - I think it's better, but I'd like to keep matching Bart's method as much as possible.

    ReplyDelete
  5. TB,
    Sorry, that's all wrong - I've just replaced the smoothed by the unsmoothed version. I'll fix that mistake. I was thinking of the other thread.

    ReplyDelete
  6. "the analysis at CA is based on the use of all CERES data...Consequently, I don't think this alternative approach is telling us anything useful about clouds."

    Please would you run the data with the re-analysis series favoured by Dessler, so we can see the difference it would make to Bart's analysis, without the Hanning filter.

    Thanks

    ReplyDelete
  7. TB,
    I can try. I mainly did the above analysis to study the method, despite feeling that the data set was wrong. Now I'm convinced the method is not useful. It returns the ratio of the trends. The sign of that depends mainly on Hadcrut3, as the CERES data has a negative trend, and it's likely that the reanalysis data will too.

    Still, it's easy enough. Not tonight, though (11pm here) and I've got a busy day tomorrow.

    ReplyDelete
  8. TB,
    I've put up a new post. The "feedback" answer came out positive, but I would not have faith in that. As I showed at CA, the sign is determined by the gradient of T. If you start at 2001 rather than 2000, it could well be negative, just as Bart's number from the CERES data changed sign.

    ReplyDelete