Recommended Authors:

MATLAB Code updated February 7, 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.

Our routines and subroutines also use an interleave function written by Jason Blackaby, and available in the public domain, but not represented here. All lines of code, with the exception of the first-line 'function' statement of each routine, as well '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)

% zero-out DC component %

M=mean(x)
B=(x)-(M)
A=B

% set N equal to size of vector A %

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

% trim end of vector A so N is odd %

if ((rem(N,2)==0))
AA=A(1:(N-1))
A=AA
N=N-1
end

% loop to trim A so (N-1) a multiple of 32 %

for g = (1:16)
NQ=((N-1)/16)
if ((rem(NQ,2)==0))
else
AA=A(1:(N-2))
A=AA
N=N-2
end
end

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

CC=forksub(A)
fork=(real(CC))+(M)
end


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

function bork = bork(x)

% flip A to achieve initial time reversal for "pre-casting" %

xf=fliplr(x)

% zero-out DC component %

M=mean(xf)
B=(xf)-(M)
A=B

% set N equal to size of vector A %

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

% trim end of A so N is odd %

if ((rem(N,2)==0))
AA=A(1:(N-1))
A=AA
N=N-1
end

% loop to trim A so (N-1) a multiple of 32 %

for g = (1:16)
NQ=((N-1)/16)
if ((rem(NQ,2)==0))
else
AA=A(1:(N-2))
A=AA
N=N-2
end
end

% call forecasting subroutine with A, undo time reversal, and add DC value back in %

CC=forksub(A)
bork=(real(fliplr(CC)))+(M)
end


Restorative Up-Sampling Master File {rup.m}
([v] being number of sampling frequency doublings, with [x] as input data {row vector})

function rup = rup(x,v)

% zero-out DC-value %

M=mean(x)
A=(x)-(M)

% 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 end of (A) so N is even %

if ((rem(N,2)==0))
else
AA=A(1:(N-1))
A=AA
N=N-1
end

% loop to ensure that for A, (N) is a multiple of 32 %

for g = (1:16)
NQ=((N)/16)
if ((rem(NQ,2)==0))
else
AA=A(1:(N-2))
A=AA
N=N-2
end

end

% call spoon subroutine with A %

B=spoon(A)
A=B
end

% take only real component and add initial DC value %

rup=real(A)+M
end


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

function forksub = forksub(x)

A=(x)

zc=complex(0.0,0.0)

Ns=size(A)

N=Ns(2:2)

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

% prepare and define extended half-length versions %

% of DFT of bandlimited input (x), giving (C) and (C2) %

Z=zeros([1 (NH-1)])

B=fft(A)

C=[B(1:NH)]

C2=[B(1:NH),Z]

% prepare and define (E), contracted double-DFT version %

% of half-length frequency signal (C2) with midsection zeroed-out %

D2=fft(C2)

D2Ns=size(D2)

D2N=D2Ns(2:2)

D2NH=(((D2N-1)/2))

D1A=D2(1:(D2NH+1))

D2A=D2((D2NH+2):D2N)

ZE=zeros([1 D2NH])

E=([D1A,ZE,ZE,D2A])

% measure length of (E) and define (ENH) %

ENs=size(E)

EN=ENs(2:2)

ENH=((EN-1)/2)

% reverse and shift to the right by (ENH), yielding (ERM) %

% then perform straight flip, giving (ERMF) %

Edc=E(1:1)

Efq=E(2:EN)

ER=[Edc,fliplr(Efq)]

ERdc=ER(1:1)

ER1=ER(2:(ENH+1))

ER2=ER((ENH+2):(EN))

ERM=[ER2,ERdc,ER1]

ERMF=fliplr(ERM)

% prepare and define (J), modified version of phantom-signal %

% which is the negative ifft of the zeroed-out midsection of (F) (below) %

% (or what it would be if filled in) as per contextual signal processing %

G=(ifft(ERMF))

GNs=size(G)

GN=GNs(2:2)

GNH=((GN-1)/2)

H=G(2:((GNH)+1))

J=[zc,H]

% multiply fft(J) by factor to give (W), and %

% use (Y1f) and (Y2f) to account for unmanifest %

% half-flip in second (modulated and somewhat in-limbo) reversal %

% yielding (X) %

W=((GN/(GNH+1))^3)*(fft(J))

WNs=size(W)

WN=WNs(2:2)

WNH=((WN-1)/2)

WNQ=((WN-1)/4)

ZW=zeros([1 WNQ])

Y1=W((1):(WNH))

Y2=W((WN-(WNH+1)+1):(WN))

Y1f=fliplr(Y1)

Y2f=fliplr(Y2)

X=[Y1f,Y2f]

% account for (undo) modulation (1) implied by difference %

% between straight-flip in (ERMF) and bonefied reversal %

XX=ifft(X)

K=(1/GN)*(linspace(GN,(GNH+1),(GNH+1)))

EX=exp(1)

L=XX.*(EX.^((1i)*(2*pi)*(-1).*(K)))

M=fft(L)

% harvest midsection components of "double-DFT %

% and order them to form Z and Zaliased %

X1=M(1:(WNQ))

X2=M((WN-WNQ+1):(WN))

X1f=fliplr(X1)

X2f=fliplr(X2)

Z=[zc,ZW,X2f,X1f,ZW]

Zaliased=[zc,X2,ZW,ZW,X1]

% prepare and define (F), contracted DFT of half-length frequency %

% signal (C) with center-section "zeroed-out" %

D=fft(C)

DNs=size(D)

DN=DNs(2:2)

DNH=(((DN-1)/2))

D1Cdc=D(1:1)

D1C=D(2:(DNH+1))

D2C=D((DNH+2):DN)

ZF=zeros([1 ((DN-1)/2)])

F=[D1Cdc,D1C,ZF,ZF,D2C]

% take (F) and add (Z) [the negative re-constructed "double-DFT" %

% of full mid-section phantom-signal, with (F)'s midsection being %

% zeroed-out, and subtract the implied aliased components of (F), %

% (Zaliased) to give fully embodied DFT of expanded %

% half-length frequency signal (M) %

M=(((DN)/(DNH))*F)+Z-Zaliased

% take ifft of (M), giving (P), and measure length %

P=ifft(M)

PNs=size(P)

PN=PNs(2:2)

% re-construct full-length frequency-domain signal, %

% take ifft to get back to forecasted time domain %

% scale and subtract DC-component %

R=([P,(fliplr(conj(P(2:PN))))])

S=ifft(R)

T=real(S)-mean(real(S))

forksub=(1/4)*(T)

end

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

function spoon = spoon(x)
zc=complex(0.0,0.0)

A=x

Ns=size(A)

N=Ns(2:2)

% section wraps (A) around unit-circle in complex plane along with reverse of (A) %

B=[zc,A,fliplr(A)]

% section takes inverse-DFT of (B) %

Cr=real(ifft(B))

% section takes half-length of ifft %

Crh=Cr(1:(N+1))

% section strips C of its mean value %

Crm=mean(Crh)

Crhn=Crh-Crm

% section forecasts Crnh to near double-length %

Er=forksub(Crhn)
ErNs=size(Er)

ErN=ErNs(2:2)

% section re-adds mean value removed earlier %

Fr=Er+Crm

% section reconstructs full-ifft %

FF=[real(Fr),real(fliplr(Fr(2:(ErN))))]

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

J=fft(FF)

% section takes first [(2*N)-1] points beginning with third point as restoratively up-sampled signal %

NU=((2*N)-1)

K=J(3:(NU+2))

% section strips away any residual DC value %

L=(real(K)-mean(real(K)))

spoon=(1/(sqrt(2)))*L

end

Quantization Error and Artifact Reduction (QEAR) algorithm
([p] being number of sampling frequency doublings, with [x] as input data {row vector})

function qear = qear(x,p)

% strip DC value %

a=x-mean(x)

% measure length of (a) %qs=size(a)
q=qs(2:2)

% set up and accomplish discrete-difference operation on (a) to yield (ad) %

a1=a(1:1)
a2=a(2:(q))
ap=[a2,a1]
ad=(ap-a)

% restoratively up-sample (ad) with (p) frequency-doublings set %

b=rup(ad,p)

% set up low-pass frequency parameters and accomplish filtration %
fsamp = (2.00^(p))
fcuts = [0.45 0.55]
mags = [1 0]
devs = [0.05 0.005]

[n,Wn,beta,ftype] = kaiserord(fcuts,mags,devs,fsamp)

hh = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale')
d=conv(b,hh)
% measure length of filter output (d) %
dns=size(d)
dn=dns(2:2)

% set up and accomplish numerical integration loop to "undo" discrete-differentiation %

pa=zeros([1 dn])
pa(1:1)=a(1:1)
for v = (2:dn)
pa(v:v)=pa((v-1):(v-1))+((trapz(d((v-1):(v)))))
end

% scale output, and downsample to original sampling frequency %

% high pass filtration eliminates low-frequency noise from numerical integration %

f=(1/(2^(p)))*downsample(pa,(2^(p)))
fr=real(f)

% high-pass filter and remove any DC value %

frhp=highpass(fr,20,22050)
qear=frhp-mean(frhp)
end