Simulated example

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.

In [1]:
%pylab inline
from GUM2DFT import GUM_DFT,GUM_iDFT,GUMdeconv, AmpPhase2DFT
from statsmodels.tsa.arima_process import ArmaProcess
Populating the interactive namespace from numpy and matplotlib

Measurement system

In [9]:
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 signal and simulated output signal

In [3]:
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)]
In [4]:
# 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))
In [8]:
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)
Out[8]:
<matplotlib.text.Text at 0x19b1d390>

Regularized inverse system for deconvolution

In [10]:
# low-pass filter
fcut = 120e6
HL = 1/(1+1j*f/fcut)**2
uAmp /= np.abs(HL)
# inverse system for input estimation    
H = FR/HL
In [12]:
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 of uncertainties

In [13]:
# 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))
In [17]:
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);
In [ ]: