Re : Signal extrapolation with FFT

[nferry]

And just a remark:

The "repeat method" is not perfect equivalent to the "fourier extrapolation method",

because this last can play with the FFT coef : (Try to decrese the numbers of coef computed)

to supress high frequency, noise filter, resolution interpolation or something else.

So the compromise between flexibility and performance is preserved. 

Nicolas

[altenbach]

---

nferry wrote:

The "repeat method" is not perfect equivalent to the "fourier extrapolation method",

because this last can play with the FFT coef

---

This thread is about plain extrapolation, so the methods are equivalent in the context of the original problem description. You did not mention filtering up until now. Filtering is entirely seperate.

Don't get me wrong, I do a lof of processing via FFTs (filtering, convolution, etc) and it is a great tool if applied to a problem that actually needs it. However, even if you use it, you never need to sum the components clumsily in a loop (unless you only want to keep very few tones). All you need to do is manipulate the complex frequency domain signal, then transform it back to get the filtered signal. Your method is exceedingly inefficient. 😉

[Chondrally]

http://ieeexplore.ieee.org/xpl/login.jsp?tp=&arnumber=651165&url=http%3A%2F%2Fieeexplore.ieee.org%2F...

Interpolation and Extrapolation using High resolution DFT by Saachi, Ulrych and Walker

[Chondrally]

In mathematica it looks something like the following.... BSD-3 open source license copyright Richard J Belshaw , ASA,IEEE,ACM,AAAS,ACS,AGU, 2014, all rights reserved

ClearAll[fourieranalysis];

fourieranalysis[actdata_List,meanpath_List,futuresteps_Integer,debugflag_Integer] := Module[{N2,N1,N3,x,coeffs,deTrendedData,volatility,x2,temp1,temp2,temp3,avg,sigman,x3,x1,x4,y,x5,y2,sigmaX,i,j,k,n,x7,x8,x9,lambda,F,A,X1,X2,X3,Q3,temp,Jcg1,count,Jcg2,X4,M,flaga},

x = actdata;

N2=2;

N1=Length[x];

N3=Length[meanpath];

While[2^N2<N1+futuresteps,N2=N2+1];

M=2^N2;

Clear[coeffs,deTrendedData,volatility];

{coeffs,deTrendedData,volatility}=getPercentErrorLinear[x[[1;;N1]]];

x2=deTrendedData;

avg = Mean[deTrendedData];

sigman =0.0;

Do[sigman=sigman+(deTrendedData[[j]]-avg)^2,{j,1,Length[deTrendedData]}];

sigman=Sqrt[sigman/(Length[deTrendedData]-1)];

x3=ConstantArray[avg,N1];

x1=ConstantArray[0.0,M];

x4=x1;

x4[[1;;Length[x2]]]=deTrendedData[[1;;Length[x2]]];

y = ForwardFFT[x4];

x5=x1;

x5[[1;;N1]]=x3;

y2 = ForwardFFT[x5];

sigmaX=0.0;

Do[sigmaX=sigmaX+(Abs[y[[i]]]-Abs[y2[[i]]])^2,{i,1,M}];

sigmaX = Sqrt[sigmaX/(M-1)];

lambda = N[sigman^2/sigmaX^2];

F = ConstantArray[0.0,{N1,M}];

Do[F[[n,k]]=N[(1/M)*Exp[Sqrt[-1]*2*Pi*(n-1)*(k-1)/N1]],{n,1,N1},{k,1,M}];

If[debugflag1,Print["*2"]];

flaga=1;

count=0;

While[flaga1&&count<3,

flaga=0;

If [N30||count>0,

x4[[N1+1;;M]]=RandomReal[{Min[x2]+.5*Abs[Mean[x2]-Min[x2]],Max[x2]-.5*Abs[Max[x2]-Mean[x2]]},M-N1],

x4[[N1+1;;N1+futuresteps]]=meanpath;

x4[[N1+futuresteps;;M]]=RandomReal[{Min[x2]+.5*Abs[Mean[x2]-Min[x2]],Max[x2]-.5*Abs[Max[x2]-Mean[x2]]},M-N1-futuresteps+1]

];

count=count+1;

If[debugflag1,Print["*4"]];

X1=ForwardFFT[x4];

If[debugflag1,Print["*5"]];

x7=x4;

If[debugflag1,Print["*5a"]];

x8 = ConstantArray[0.0,M];

x8[[1;;N1]]=x2;

x8[[N1+1;;M]]=x4[[N1+1;;M]]*1.01;

X2=ForwardFFT[x8];

If[debugflag1,Print["*9"]];

temp1 = x8[[N1+1;;M]];

temp2 = x7[[N1+1;;M]];

Jcg2=N[Abs[1-temp1/temp2]];

temp = matlabFind1[Jcg2,Greater,10^(-4)];

j=1;

Q3=ConstantArray[0.0,{M,M}];

x9=ConstantArray[0.0,M];

While[Length[temp]>0&&j<M-N1,

X3=ForwardFFT[x8];

Do[Q3[[i,i]]=1+Conjugate[X3[[i]]]*X3[[i]]/sigmaX^2,{i,1,M}];

A=lambda*IdentityMatrix[N1]+N[F.Q3.ConjugateTranspose[F]];

X4= N[Q3.ConjugateTranspose[F].Inverse[A].x2];

(*If [j1,

maxval = Max[Abs[X4[[1;;M/2]]]];

posn = matlabFind1[Abs[X4[[2;;M/2]]],GreaterEqual,maxval];

X4[[posn+1]]=2*X4[[posn+1]];

X4[[M-posn+1]]=2*X4[[M-posn+1]];];*)

x8=Re[IFFT[X4]];

x8[[1;;N1]]=x2;

X4=ForwardFFT[x8];

Jcg1=Jcg2;

temp1 = x8[[N1+1;;M]];

temp2 = x7[[N1+1;;M]];

Jcg2=N[Abs[1-temp1/temp2]];

temp = matlabFind1[Jcg2,Greater,10^(-4)];

If [Length[temp]>0,

temp3= temp1[[temp]]-Jcg2[[temp]]*(temp1[[temp]]-temp2[[temp]])/(Jcg2[[temp]]-Jcg1[[temp]]);

x9 = x8;

x9[[N1+temp]]=temp3;

x7=x8;

x8=x9;

If[debugflag1,Print["j= ",j]];

j=j+1;

X1=X2;

X2=X3;

X3=X4;

];

];

If[jM-N1,flaga=1];

];

If[debugflag1,Print["*1"]];

x9[[1;;N1]]=deTrendedData[[1;;N1]];

If[debugflag1,Print["*2"]];

x9[[N1+1;;M]]=2*x9[[N1+1;;M]];

If[debugflag1,Print["*3"]];

x8=x9[[1;;M]];

If[debugflag1,Print["*4"]];

Do[x8[[j]]=x8[[j]]+coeffs[[1]]*j+coeffs[[2]],{j,1,M}];

{x8,x8[[N1+1;;N1+futuresteps]]}]