Recommended Authors:

MATLAB Code - Update - May 24, 2023

All of the following MATLAB files must be available to the MATLAB program in order for them to work together in concert. There are a number of interdependent calls of routines and subroutines in running any of the master Fourier forecasting or restorative up-sampling files, so they all must be installed together. Admittedly, our documentation is somewhat sparse.

The {forksub.m} routine makes use of an interleave function, written by Jason Blackaby. This code is not included on this site, but it is available in the public domain. All lines of code, with the exception of the first-line 'function' statement of each routine, as well comments and 'if,' 'else,' and 'for' statements and their corresponding 'end' statements, should have semi-colons placed at their endings, a punctuation mark which this web-builder does not support.

Forwards Forecasting Master File {fork.m}
([x] as input data {row vector})

function fork = fork(x)

% declare (A) %

A=(x)-mean(x)

% set N equal to size of vector A %

Ns=size(A)
N=Ns(2:2)

% trim end of vector A so N is even %

if ((rem(N,2)==0))
else

AA=A(1:(N-1))
A=AA
N=N-1
end

% call forecasting subroutine with A, then add DC value back in %

AA=forksub(A)

BB=AA(1:N)

CC=AA*(rms(A)/rms(BB))

fork=CC

end


Backwards Forecasting Master File {bork.m}
([x] as input data {row vector})

function bork = bork(x)

% declare (A) %

A=(x)-mean(x)

% set N equal to size of vector A %

Ns=size(A)
N=Ns(2:2)

% trim end of A so N is even %

if ((rem(N,2)==0))
else

AA=A(1:(N-1))
A=AA
N=N-1
end

% call backwards-forecasting subroutine with A, and scale %

AA=borksub(A)

BB=AA((N+1):(2*N))

CC=AA*(rms(A)/rms(BB))

bork=CC
end

Restorative Up-Sampling Master File {rup.m}
([v] being number of sampling frequency doublings, [f] being the number of forecasting "passes"

and with [x] as input data {row vector})

function rup = rup(x,v,f)

% zero-out DC-value %

M=mean(x)

A=(x)-(M)

AAA=A

% loop for v frequency doublings %

for p = (1:v)

% set N equal to length of vector (A) %

Ns=size(A)

N=Ns(2:2)

% trim A so (N) is odd %

if ((rem((N-1),2)==0))

else

AA=A(1:(N-1))

A=AA

N=N-1

end

% if (f=0), define (FS), otherwise, initialize (F) and (B) for forecasting %

if (f==0)

FS=A

else

F=A(1:(N-1))

B=A(2:(N))

end

% loop for accomplishes (f) number of pre- and post- forecasts %

if (f==0)

else

for rr=(1:f)

NFs=size(F)

NF=NFs(2:2)

FF=fork(F))

BB=bork(B)

F=F

B=B

end

% prepare signal - sandwiched between pre- and post- forecasts %

BNs=size(B)

BN=BNs(2:2)

BS=B(1:(BN-(N)))

FS=[BS,F]

F=FS

B=FS

end

% call spoon subroutine with (FS) and then trim to "area of concern" %

BU=spoon(FS)

BUNs=size(BU)

BUN=BUNs(2:2)

BUNH=(((BUN-1)/2)+1)

FUN=BU((BUNH-(N-1)):(BUNH+(N-1)))

A=FUN

end

% take only real component and scale %

rup=(real(A))*(rms(AAA)/rms(A))

end

Forward Forecasting Subroutine File (forksub.m) (subroutine)

function forksub = forksub(x)

% initialize some variables %

zc=complex(0.0,0.0)

A=fft(x)

B=(x)

% measure length (BN) of input, (B) %

% and define (BNH) and zero array (ZB) %

BNs=size(B)

BN=BNs(2:2)

BNH=(BN/2)

ZB=zeros([1 (BN-1)])

% accomplish conventional (x4) up-sampling of (B) by zero-padding %

% ib the frequency-domain to give (B4) %

BF=fft(B)

BF1=BF(1:BNH)

BF2=BF((BNH+2):BN)

BFM=BF((BNH+1):(BNH+1))

BFU=2*([BF1,(0.5)*BFM,ZB,(0.5)*BFM,BF2])

C=real(ifft(BFU))

% measure length of (C) %

CNs=size(C)

CN=CNs(2:2)

% add "nulls" to (C), yielding (D) %

Z1=zeros([1 (CN)])

D=interleave(Z1,C)

% take DFT of (D) %

F=fft(D)

% expand (F) by interleaving with zeros %

Z2=zeros([1 (2*CN)])

Gr=(interleave((real(F)),Z2))

Gi=(interleave((imag(F)),Z2))

G=((complex((Gr'),(Gi')))')

% return to time-domain and negate "null-values" %

H=-ifft(G)

% measure length of (H), and declare half-length (HNH) %

HNs=size(H)

HN=HNs(2:2)

HNH=(HN/2)

% take only first half of (H), yielding (J) %

J=H(1:HNH)

% down-sample (J) by factor of two, with one-point offset, yielding (JJ) %

JJ=downsample(J,2,1)

% take DFT of (JJ), and down-sample halves of (M) %

% to give odd-frequency components (MM) %

M=fft(JJ)

MNs=size(M)

MN=MNs(2:2)

MNH=((MN)/2)

M1=(M(1:(MNH+1)))

M2=(M((MNH+2):(MN)))

MM1=downsample(M1,2,1)

MM2=downsample(M2,2)

MM=[MM1,MM2]

% interleave ("to be" even) [original (A)] frequency components %

% and derived odd-frequency components (MM), take inverse DFT and %

% scale to give non-bandlimited forecast %

Pr=interleave(real(A),real(MM))

Pi=interleave(imag(A),imag(MM))

P=((complex((Pr'),(Pi')))')

forksub=(2)*(real(ifft(P)))

end

Backward Forecasting Subroutine File (borksub.m) (subroutine)

function borksub = borksub(x)

% initialize some variables %

Ns=size(x)

N=Ns(2:2)

B=[(x(1:1)),fliplr(x(2:N))]

A=fft(B)

% measure length (BN) of input, (B) %

% and define (BNH) and zero array (ZB) %

BNs=size(B)

BN=BNs(2:2)

BNH=(BN/2)

ZB=zeros([1 (BN-1)])

BF=fft(B)

BF1=BF(1:BNH)

BF2=BF((BNH+2):BN)

BFM=BF((BNH+1):(BNH+1))

BFU=4*([BF1,(0.5)*BFM,ZB,(0.5)*BFM,BF2])

C=real(ifft(BFU))

% measure length of (C) %

CNs=size(C)

CN=CNs(2:2)

% add even-numbered "nulls" to (C) %

Z1=zeros([1 (CN)])

D=interleave(Z1,C)

% take DFT of (D) %

F=fft(D)

% expand (F) by interleaving with zeros %

Z2=zeros([1 (2*CN)])

Gr=(interleave((real(F)),Z2))

Gi=(interleave((imag(F)),Z2))

G=((complex((Gr'),(Gi')))')

% return to time-domain and negate "null-values" %

H=-ifft(G)

% measure length of (H), and declare half-length (HNH) %

HNs=size(H)

HN=HNs(2:2)

HNH=(HN/2)

% take only first half of (H), yielding (J) %

J=H(1:HNH)

% down-sample (J) by factor of two, with one-point offset, yielding (JJ) %

JJ=downsample(J,2,1)

% take DFT of (JJ), and down-sample halves of (M) %

% to give odd-frequecy components (MM) %

M=fft(JJ)

MNs=size(M)

MN=MNs(2:2)

MNH=((MN)/2)

M1=(M(1:(MNH+1)))

M2=(M((MNH+2):(MN)))

MM1=downsample(M1,2,1)

MM2=downsample(M2,2)

MM=[MM1,MM2]

% interleave ("to be" even) [original (A)] frequency components %

% and derived odd-frequency components (MM), take inverse DFT and %

% scale to give non-bandlimited forecast %

Pr=interleave(real(A),real(MM))

Pi=interleave(imag(A),imag(MM))

P=((complex((Pr'),(Pi')))')

Z=(2)*(real(ifft(P)))

borksub=[fliplr(Z(2:(2*N))),Z(1:1)]

end

Up-Sampling Subroutine File (spoon.m) (subroutine)

function spoon = spoon(x)

% initialize some variables %

Ns=size(x)

N=Ns(2:2)

NH=(((N-1)/2)+1)

% section arranges input as if it were in the frequency-domain %

A=real([(x(NH:N)),fliplr(x(1:(NH-1)))])

% section takes inverse-DFT of (A) %

Cr=real(ifft(A))

Ci=imag(ifft(A))

% section ensures portion to be forecasted is even length %

if ((rem(NH,2))==0)

Dr=Cr(1:(NH))

Di=Ci(1:(NH))

else

Dr=Cr(1:(NH-1))

Di=Ci(1:(NH-1))

end

% section strips (Dr) and (Di) of their mean value %

Drm=mean(Dr)

Drn=Dr-Drm

Dim=mean(Di)

Din=Di-Dim

% section forecasts (Dim) and (Din) twice %

Er=forksub(forksub(Drn))

Ei=forksub(forksub(Din))

% section trims forcast to needed length %

Fr=Er(1:((2*NH)-1))

Fi=Ei(1:((2*NH)-1))

% section reconstructs full-length inverse DFT %

Gr=[(real(Fr+(0.5*Drm))),(real(fliplr(Fr(2:((2*NH)-1))+(0.5*Drm))))]

Gi=[(real(Fi+(0.5*Dim))),(-real(fliplr(Fi(2:((2*NH)-1))+(0.5*Dim))))]

G=Gr+((1i)*Gi)

% section returns to "time-domain" by taking a DFT %

J=fft(G)

% section reconstructs [(2*N)-1] point time-domain up-sampled signal %

L=real([fliplr(J((N+1):((2*N)-1))),(J(1:N))])

% section strips away any residual DC value and declares final variable %

M=(L)-mean(L)

spoon=M

end

Quantization Error and Artifact Reduction (QEAR) algorithm

function qear = qear(x,p,v,f)

% initialize (A) %

A=(x)

% loop for number of passes (p) %

for z=(1:p)

% measure length of (A) %

Ns=size(A)

N=Ns(2:2)

% trim (A) to odd-numbered length %

if ((rem(A,2))==0)

AA=A(1:(N-1))

A=AA

N=N-1

end

% set up and accomplish discrete-difference for (A) %

AS=[A(N:N),A(1:(N-1))]

D=A-AS

% restoratively up-sample (D) with %

% (v)-frequency doublings and (f) forecast loops %

E=rup(D,v,f)

% perform simple low-pass filtration of (E) %

F=lowpass(E,0.5,(2^v))

% measure length of (F) and set-up %

% (G) to receive values of numerical integration loop %

FNs=size(F)

FN=FNs(2:2)

G=zeros([1 FN])

% accomplish numerical integration %

for y=(2:FN)

G(y:y)=G((y-1):(y-1))+F(y:y)

end

% down-sample (G) back to original sampling frequency %

H=downsample(G,(2^v))

% perform simple highpass filtration %

% to remove low-frequency integration noise %

J=highpass(H,0.0001,1)

end

% scale and define output (qear) of noise reduction algorithm %

Z=J*(rms(A)/rms(J))

qear=real(Z)

end