This is a clean simulated deconvolution example. That is, using the measured frequency response an output signal is generated by means of convolution and thereafter the input signal estimated by deconvolution. The advantage of this approach is to verify the general applicability of the method.
%pylab inline
from GUM2DFT import GUM_DFT,GUM_iDFT,GUMdeconv, AmpPhase2DFT
from statsmodels.tsa.arima_process import ArmaProcess
calib = np.loadtxt("data/calib_PA055.dat")
f = calib[:,0]
Nf = 2*len(f)-1
FR = calib[:,1]*np.exp(1j*calib[:,3])
uAmp = calib[:,2]
uPhas= calib[:,4]
input_file = "data/reference_signal.dat"
noise_std = 1e-4
time,pressure = np.loadtxt(input_file).T
Ts = 2e-9
P = np.fft.rfft(np.r_[pressure,np.zeros(2*len(f)-len(pressure)-1)])
x = np.fft.irfft(FR*P)[:len(time)]
# adding correlated noise (ARMA process)
N = len(time)
ar = np.array([0.75,-0.6])
ma = np.array([0.45,0.1])
ARMA = ArmaProcess(ar,ma)
noise = ARMA.generate_sample(N)*0.01
x += np.random.randn(len(time))*noise_std
acov = ARMA.acovf(nobs=N)*noise_std**2
Ux = np.zeros((N,N))
for k in range(N):
Ux[k,k:] = acov[:N-k]
Ux = Ux + Ux.T - np.diag(np.diag(Ux))
figure(figsize=(14,8));clf()
plot(time*1e6,pressure,label="true pressure signal")
plot(time*1e6,x,label="simulated measurement")
legend(loc="best");
xlabel(r"time / $\mu$s",fontsize=20)
ylabel("pressure / MPa",fontsize=20)
# low-pass filter
fcut = 120e6
HL = 1/(1+1j*f/fcut)**2
uAmp /= np.abs(HL)
# inverse system for input estimation
H = FR/HL
figure(figsize=(14,12));clf()
subplot(211)
errorbar(f*1e-6,calib[:,1],2*uAmp,fmt=".",alpha=0.2)
xlim(0.5,150);ylim(0,0.25)
xlabel("frequency / MHz",fontsize=22); tick_params(which="both",labelsize=18)
ylabel("amplitude / a.u.",fontsize=22)
subplot(212)
errorbar(f*1e-6,calib[:,2],2*uPhas,fmt=".",alpha=0.2)
xlim(0.5,150)
xlabel("frequency / MHz",fontsize=22); tick_params(which="both",labelsize=18)
ylabel("phase / deg",fontsize=22);
# propagation to real/imag domain
UH = AmpPhase2DFT(np.abs(FR),np.angle(FR)*np.pi/180,np.r_[uAmp,uPhas*np.pi/180])[1]
# propagation from time to frequency domain
Y,UY = GUM_DFT(x,Ux,N=Nf,return_complex=True)
# propagation through deconvolution equation
XH,UXH = GUMdeconv(H,Y,UH,UY)
# propagation back to time domain
xh,Uxh = GUM_iDFT(XH,UXH,Nx=len(x))
ux = np.sqrt(np.diag(Uxh))
figure(figsize=(14,8));clf()
plot(time*1e6,xh,label="estimate of measurand",linewidth=2)
plot(time*1e6,pressure,"--",label="true value of measurand",linewidth=2)
fill_between(time*1e6,xh+ux,xh-ux,alpha=0.2)
xlabel(r"time / $\mu$s",fontsize=22)
ylabel("signal amplitude / a.u.",fontsize=22)
tick_params(which="major",labelsize=18)
legend(loc="best",fontsize=20,fancybox=True)
xlim(0.5,1.5);