# How to construct optimal Wiener Filter?

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.

• In your third plot at x=0, why is P_S + P_N (green, -3) less than P_S alone (red, +5)? And why does one of the five things in the legend box have "ln()"? Is there a chance you're mixing logarithmed and unlogarithmed things together in your algebra (additions, subtractions)? While the log of P_N can be negative, P_N will be small but always positive, so P_S + P_N should always be larger than P_S alone.
– uhoh
Commented May 16, 2023 at 10:34
• In fact, isn't there something in the old testament of the bible about "mixing together of logarithmed and unlogarithmed things"? If not, it seems like there should be (humor!)
– uhoh
Commented May 16, 2023 at 10:38
• you pointed out the mistake, I mixed up 2 variable names in python. I edited the question with the correct results Commented May 16, 2023 at 11:39
• okay see my comment under my answer; I've rolled back the edit (the text is still available) so you can paste it into a new answer post. Thanks!
– uhoh
Commented May 16, 2023 at 12:09

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()

• I work in space of natural logarithm but it should not make a difference for the approach. I added a plot where I included the signal and noise models extrapolated over the whole spectrum, as well as their sum. The problem is that, the sum does not reproduce the actual powerspectrum (because the function values are negative). Commented May 15, 2023 at 8:58
• Also I am not sure about : " It is the complete function W(k) that you least-squares fit to the signal over the full range of k". I fitted the signal and model part for the corresponding regions using scipy curve_fit()', so it extends over the full range of k if I bring the 2 parts together. Is that what was meant in the quote ? ! Commented May 15, 2023 at 9:00
• @trynerror hmm... well I did my best to answer the question as it was written. I mean that you should only do one fit, including everything in it, not in two parts. But I suppose you could do two parts separately if you want. Now that the question has been changed I can't tell if my answer is still helpful. Question marks appear in your title, and in the body of your question as "Does anyone have an idea how to apply this filter equation to this case?" Have I answered these questions?
– uhoh
Commented May 15, 2023 at 10:42
• Yes, your answer was helpful, I just edited the question to point out the problem I still have. See the newest edit, my green curve in the 3rd plot does not reproduce the form of the powerspectrum very well, and I am trying to figure out how to improve that Commented May 15, 2023 at 10:45
• @trynerror No. That's not how Stack Exchange works. We do not keep editing a question post to add the solution. If you have resolved the issue, the only place you should consider putting that is in a comment or better yet an new answer post. I've rolled back your edit. The text is still there (just click "edited" under the question) and you can copy/paste that to new answer post. thanks!
– uhoh
Commented May 16, 2023 at 12:04

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.

• Yay, congrats on working the problem and finding the solution!
– uhoh
Commented May 16, 2023 at 20:53