How to construct optimal Wiener Filter?

Question

I am trying to construct a Wiener Filter, to filter the ratio of the peak from the cross-correlation function, between a galaxy spectra and a template spectra, with the peak of the auto-correlation function of the template. Therefor, I calculated the powerspectrum P(k) of the cross-correlation peak in Fourier-Space and normalized it (by deviding through its maximum): $P(k)=|X_{Peak}(k)|^2$. This powerspectrum decomposes into a noise part and a signal part which become clearly visible in the logarithmic space. The following plot shows $\ln{P(k)}$:

In eq. (32) of Brault & White (1971) (https://ui.adsabs.harvard.edu/abs/1971A%26A....13..169B/abstract) the Wiener Filter is described as: $W(k) = \frac{P'_S(k)}{P'_S(k)+P'_N(k)}$, where $P'_S(k)$ is a simple model for the signal part and $P'_N(k)$ a simple model for the noise part. I chose a parabola for the signal (central) part and linear functions (left and right tails) for the noise. With the following result:

Since I work with numeric arrays I am not sure how to actually apply eq. (32) for the Wiener Filter, since $P'_S$ is of course a shorter array than $P'_S(k)+P'_N(k)$, so the numerical division fails. Does anyone have an idea how to apply this filter equation to this case? After constructing the Wiener-Filter in the logarithmic space, I wanted to transfer it back using $\exp{W(k)}$.

When I extrapolate my signal and noise model(and fit only in the appropriate reason) to the full range it looks like this. The green curve shows P_S+P_N, where I calculated the sum in non-logarithmic space and then transferred back using np.log(). This green curve does not look like the sum in the given answer. The orange one is still the model for the spectrum ( concatenating the arrays). When I just follow the steps and calculate W_log = P_S /(P_S+P_N) the problem is that all my function values are negative, that is why I first use np.exp() before adding them together. The resulting Wiener filter I get, does not look right, I also plotted the ratio of correlation peaks which this Filter will get multiplied by.

Answer 1

update: This may not yet be the complete answer, stay tuned...

---

tl;dr:

To make this work and get the smooth "knee" near the "cutoff" at $k = \pm 0.005$ employ the oft-used trick of adding two functions with very different behavior such that each grows much larger than the other moving away from the cutoff in opposite directions and putting it in the denominator. Thus each of your functions must be evaluated over the full range of $k$:

---

Looking at Fig. 15 in your reference, it seems to be that even though your parabolic signal and linear noise models dominate in different regions, for both the fitting step and the Wiener filter generating step you should extend each part of the function arrays to the full range of $-0.017 \le k \le +0.017$ (or $0 \le k \le +0.017$) and then perform the division.

For example, in your plot it looks like you are using something like

$$ \log_{10}(P'_S(k)) = -a k^2 $$ $$\log_{10}(P'_N(k)) = -b |k| - c$$

or

$$P'_S(k) = 10^{-a k^2} $$ $$P'_N(k) = 10^{-b |k| - c}$$

with (very roughly) the following coefficients

$$a, b, c = -6 \times 10^5, -800, -11$$

Don't forget to take the absolute value $|k|$ for the $b$ term of the noise!

Then:

$$W(k) = \frac{P'_S(k)}{P'_S(k) + P'_N(k)}$$

NOTE: It is the complete function $W(k)$ that you least-squares fit to the signal over the full range of $k$!

(If you need help with that part let me know!)

The Python script below suggests how to manage logarithms, division etc.

Fig. 15 Illustration of (A) how the smooth power models of the signal and noise are determined by the least squares fitting of the computed power spectrum of the data, and then used to construct the optimum filter as shown in (B)

---

Generating script:

import numpy as np
import matplotlib.pyplot as plt

k = np.linspace(-0.017, 0.017, 3401)

a, b, c = -6E+05, -800, -11

log10_Ps = a * k**2
log10_Pn = b * np.abs(k) + c

Ps = 10**log10_Ps
Pn = 10**log10_Pn
P_sum = Ps + Pn

W = Ps / P_sum

thingz = (Ps, Pn, P_sum), (W, )
namez = ('Ps', 'Pn', 'Ps+Pn'), ('W=Ps/(Ps+Pn)', )

fig, axes = plt.subplots(2, 1)
ax1, ax2 = axes

for ax, things, names in zip(axes, thingz, namez):
    for thing, name in zip(things, names):
        ax.plot(k, thing, label=name)
    ax.set_yscale('log')
    ax.set_ylim(1E-25, 1E+2)
    ax.legend(facecolor='none', edgecolor='none', fontsize=14)
ib = np.argmax(k >= 0.005) # box
for thing, ax in zip((P_sum, W), axes):
    ax.plot(k[ib:ib+1], thing[ib:ib+1], 'sr',
            markersize=16, fillstyle='none')
# axes[1].set_xlim(0, 0.01)
plt.show()

if True:
    fig, ax = plt.subplots(1, 1)
    for thing, name in zip(thingz[0], namez[0]):
        ax.plot(k, thing, label=name)
    ax.set_yscale('log')
    ax.set_xlim(0.0047, 0.0055)
    ax.set_ylim(1E-16, 1E-13)
    ax.legend(facecolor='none', edgecolor='none', fontsize=14)
    plt.show()

Answer 2

To construct the optimal Wiener Filter $W(k)$:

1.) compute power spectrum of the crosscorrelation(between a galaxy spectrum and a stellar template spectrum peak) $X_{G, S}(z) = \int dx \, G(x+z)S(x)$, where $G(x)$ is the galaxy, and $S(x)$ is the template. The peak should be ca. 4 times the FWHM of $X_{G, S}$. The powerspectrum of the crosscorrelation peak is: $P= |\tilde{X}_{G, S, Peak}|^2$, where $\tilde{X}_{G, S, Peak}$ is the fourier transform of $X_{G, S, Peak}$.

Python code to obtain $X_{G, S}$:

corr_gal_tem = correlate(flux_gal, flux_temp, mode='same')

cut out the peak to obain $X_{G, S, Peak}$:

corr_gal_tem_peak = corr_gal_tem[peak_index_ccf-cut_off_pixel_2-1:peak_index_ccf+cut_off_pixel_2]

taper ends of functions for fourier transform:

corr_peak_taper = np.kaiser(len(corr_gal_tem_peak), 8)
corr_gal_tem_peak = corr_gal_tem_peak * corr_peak_taper

perform fourier transform to get $\tilde{X}_{G, S, Peak}$:

corr_gal_tem_peak_fourier = np.fft.fftshift(fourier(corr_gal_tem_peak))

and the power spectrum $P$:

powspec_gal_tem = np.power(np.abs(corr_gal_tem_peak_fourier), 2)

2.) The Wiener Filter is defined by Brault and White (1971): $W(k)=\frac{P'_S(k)}{P'_S(k)+P'_N(k)}$, where $P'_S(k)$ is a model for the signal part and $P'_N(k)$ a simple model for the noise part. In logarithmic space, the signal part and the noise part can be modeled by a parabola and linear functions respectively. In the following plot you can see the fitted and extrapolated functions in logarithmic space of the powerspectrum:

The fit is only applied in the region of the signal or noise part, therefore one has to define a cut of frequency between signal and noise. After that, one extrapolates the fitted models over the whole frequency range. These same length arrays are transferred back to non-logarithmic space and added there! After that one can apply the given Wiener formula.

    s_inter = 0.0028 # cutoff frequency

    # signal part
    s_signal = s[(s >= -s_inter) & (s <= s_inter)]
    s_signal_left = s_signal[(s_signal <=0)]
    s_signal_right = s_signal[(s_signal >0)]
    log_powspec_signal = np.log(powspec_gal_tem)[(s >= -s_inter) & (s <= s_inter)]

    # noise part
    s_noise_left = s[(s < -s_inter)]
    log_powspec_noise_left = np.log(powspec_gal_tem)[(s < -s_inter)]
    s_noise_right = s[(s > s_inter)]
    log_powspec_noise_right = np.log(powspec_gal_tem)[(s > s_inter)]
    signal_model_params = curve_fit(signal_model, s_signal, 
    log_powspec_signal)[0]
    signal_model_vals = signal_model(s_signal, *signal_model_params)
    noise_model_left_params = curve_fit(noise_model, s_noise_left, 
    log_powspec_noise_left)[0]
    noise_model_left_vals = noise_model(s_noise_left, 
    *noise_model_left_params)
    noise_model_right_params = curve_fit(noise_model, s_noise_right, 
    log_powspec_noise_right)[0]
    noise_model_right_vals = noise_model(s_noise_right, 
   *noise_model_right_params)

    #extrapolate values of signal to left and right parts
    signal_model_left = signal_model(s_noise_left, *signal_model_params)
    signal_model_right = signal_model(s_noise_right, *signal_model_params)
    signal_model_all = np.concatenate((signal_model_left, signal_model_vals, 
    signal_model_right), axis=0)
    #extrapolate values of the noise to the center
    noise_model_center_left = noise_model(s_signal_left, 
    *noise_model_left_params)
    noise_model_center_right = noise_model(s_signal_right, 
    *noise_model_right_params)
    noise_model_all = np.concatenate((noise_model_left_vals, 
    noise_model_center_left, noise_model_center_right, 
    noise_model_right_vals), axis=0)

    # add signal and noise in non logarithmic space
    P_signal = np.exp (signal_model_all)
    P_noise = np.exp(noise_model_all)
    P_sum = P_signal+P_noise

    # construct wiener filter: W(k) = P_S(k) /( P_S(k)+P_N(k))
    wiener = P_signal/P_sum

There are probably many ways how you can implement this shorter. The Wiener Filter finally looks like this, the other function is the ratio of correlation peaks and to be filtered.