<p><strong>update:</strong> This may not yet be the complete answer, stay tuned...</p>
<hr />
<h3>tl;dr:</h3>
<p>To make this work and get the smooth &quot;knee&quot; near the &quot;cutoff&quot; at <span class="math-container">$k = \pm 0.005$</span> 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. <em><strong>Thus each of your functions must be evaluated over the full range of <span class="math-container">$k$</span>:</strong></em></p>
<p><a href="https://i.sstatic.net/Wrp5b.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/Wrp5b.png" alt="plot of Wiener filter implementation from the python script below" /></a></p>
<hr />
<p>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 <strong>both</strong> the fitting step and the Wiener filter generating step you should extend each part of the function arrays to the full range of <span class="math-container">$-0.017 \le k \le +0.017$</span> (or <span class="math-container">$0 \le k \le +0.017$</span>) <em>and then</em> perform the division.</p>
<p>For example, in your plot it looks like you are using something like</p>
<p><span class="math-container">$$ \log_{10}(P'_S(k)) = -a k^2 $$</span>
<span class="math-container">$$\log_{10}(P'_N(k)) = -b |k| - c$$</span></p>
<p>or</p>
<p><span class="math-container">$$P'_S(k) = 10^{-a k^2} $$</span>
<span class="math-container">$$P'_N(k) = 10^{-b |k| - c}$$</span></p>
<p>with (very roughly) the following coefficients</p>
<p><span class="math-container">$$a, b, c = -6 \times 10^5, -800, -11$$</span></p>
<p><em><strong>Don't forget</strong></em> to take the absolute value <span class="math-container">$|k|$</span> for the <span class="math-container">$b$</span> term of the noise!</p>
<p>Then:</p>
<p><span class="math-container">$$W(k) = \frac{P'_S(k)}{P'_S(k) + P'_N(k)}$$</span></p>
<p><strong>NOTE: It is the complete function <span class="math-container">$W(k)$</span> that you least-squares fit to the signal over the full range of <span class="math-container">$k$</span>!</strong></p>
<p>(<em>If you need help with that part let me know!</em>)</p>
<p>The Python script below suggests how to manage logarithms, division etc.</p>
<p><a href="https://i.sstatic.net/Gv51H.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/Gv51H.png" alt="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)" /></a></p>
<blockquote>
<p>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)</p>
</blockquote>
<hr />
<p><a href="https://i.sstatic.net/wz7hT.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/wz7hT.png" alt="plot of Wiener filter implementation from the python script below" /></a></p>
<p>Generating script:</p>
<pre><code>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 &gt;= 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()
</code></pre>

CC BY-SA 4.0
2023-05-15T08:58:31.263Z
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).
trynerror
CC BY-SA 4.0
2023-05-15T09:00:40.380Z
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 ?
!
trynerror
CC BY-SA 4.0
2023-05-15T10:42:16.227Z
@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
CC BY-SA 4.0
2023-05-15T10:45:49.630Z
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
trynerror
CC BY-SA 4.0
2023-05-16T10:41:09.783Z
@trynerror thanks for the "accept" vote but please consider unaccepting for now until we get this all figured out. It may turn out that you find a bug yourself, in which case it's perfectly okay to post an additional answer to to your own question, and even a few days later accepting it. Goals is always to get things solved :-)
uhoh
CC BY-SA 4.0
2023-05-16T11:38:01.197Z
I will edit my question in the evening, to show how I resolved the problem.
trynerror
CC BY-SA 4.0
2023-05-16T12:04:14.977Z
@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`](https://astronomy.stackexchange.com/posts/53681/revisions)" under the question) and you can copy/paste that to new answer post. *thanks!*
uhoh
CC BY-SA 4.0
2023-05-14T22:35:16.740Z
2023-05-16T10:42:03.430Z
2023-05-16T10:42:03.430Z
uhoh
<p>To construct the optimal Wiener Filter <span class="math-container">$W(k)$</span>:</p>
<p>1.) compute power spectrum of the crosscorrelation(between a galaxy spectrum and a stellar template spectrum peak) <span class="math-container">$X_{G, S}(z) = \int dx \, G(x+z)S(x)$</span>, where <span class="math-container">$G(x)$</span> is the galaxy, and <span class="math-container">$S(x)$</span> is the template. The peak should be ca. 4 times the FWHM of <span class="math-container">$X_{G, S}$</span>. The powerspectrum of the crosscorrelation peak is: <span class="math-container">$P= |\tilde{X}_{G, S, Peak}|^2$</span>, where <span class="math-container">$\tilde{X}_{G, S, Peak}$</span> is the fourier transform of <span class="math-container">$X_{G, S, Peak}$</span>.</p>
<p>Python code to obtain <span class="math-container">$X_{G, S}$</span>:</p>
<pre><code>corr_gal_tem = correlate(flux_gal, flux_temp, mode='same')
</code></pre>
<p>cut out the peak to obain <span class="math-container">$X_{G, S, Peak}$</span>:</p>
<pre><code>corr_gal_tem_peak = corr_gal_tem[peak_index_ccf-cut_off_pixel_2-1:peak_index_ccf+cut_off_pixel_2]
</code></pre>
<p>taper ends of functions for fourier transform:</p>
<pre><code>corr_peak_taper = np.kaiser(len(corr_gal_tem_peak), 8)
corr_gal_tem_peak = corr_gal_tem_peak * corr_peak_taper
</code></pre>
<p>perform fourier transform to get <span class="math-container">$\tilde{X}_{G, S, Peak}$</span>:</p>
<pre><code>corr_gal_tem_peak_fourier = np.fft.fftshift(fourier(corr_gal_tem_peak))
</code></pre>
<p>and the power spectrum <span class="math-container">$P$</span>:</p>
<pre><code>powspec_gal_tem = np.power(np.abs(corr_gal_tem_peak_fourier), 2)
</code></pre>
<p>2.) The Wiener Filter is defined by Brault and White (1971):
<span class="math-container">$W(k)=\frac{P'_S(k)}{P'_S(k)+P'_N(k)}$</span>, where <span class="math-container">$P'_S(k)$</span> is a model for the signal part and <span class="math-container">$P'_N(k)$</span> 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:</p>
<p><a href="https://i.sstatic.net/SBwup.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/SBwup.png" alt="enter image description here" /></a></p>
<p>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.</p>
<pre><code>    s_inter = 0.0028 # cutoff frequency

    # signal part
    s_signal = s[(s &gt;= -s_inter) &amp; (s &lt;= s_inter)]
    s_signal_left = s_signal[(s_signal &lt;=0)]
    s_signal_right = s_signal[(s_signal &gt;0)]
    log_powspec_signal = np.log(powspec_gal_tem)[(s &gt;= -s_inter) &amp; (s &lt;= s_inter)]

    # noise part
    s_noise_left = s[(s &lt; -s_inter)]
    log_powspec_noise_left = np.log(powspec_gal_tem)[(s &lt; -s_inter)]
    s_noise_right = s[(s &gt; s_inter)]
    log_powspec_noise_right = np.log(powspec_gal_tem)[(s &gt; 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

</code></pre>
<p>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.
<a href="https://i.sstatic.net/MGXqn.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/MGXqn.png" alt="enter image description here" /></a></p>

CC BY-SA 4.0
2023-05-16T20:53:25.223Z
Yay, congrats on working the problem and finding the solution!
uhoh
CC BY-SA 4.0
2023-05-16T20:30:53.933Z
2023-05-16T20:30:53.933Z
trynerror
<p>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):
<span class="math-container">$P(k)=|X_{Peak}(k)|^2$</span>. This powerspectrum decomposes into a noise part and a signal part which become clearly visible in the logarithmic space. The following plot shows <span class="math-container">$\ln{P(k)}$</span>:
<a href="https://i.sstatic.net/zHuhh.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/zHuhh.png" alt="enter image description here" /></a></p>
<p>In eq. (32) of Brault &amp; White (1971) (<a href="https://ui.adsabs.harvard.edu/abs/1971A%26A....13..169B/abstract" rel="nofollow noreferrer">https://ui.adsabs.harvard.edu/abs/1971A%26A....13..169B/abstract</a>) the Wiener Filter is described as: <span class="math-container">$W(k) = \frac{P'_S(k)}{P'_S(k)+P'_N(k)}$</span>, where <span class="math-container">$P'_S(k)$</span> is a simple model for the signal part and <span class="math-container">$P'_N(k)$</span> 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:</p>
<p><a href="https://i.sstatic.net/zAKAr.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/zAKAr.png" alt="enter image description here" /></a></p>
<p>Since I work with numeric arrays I am not sure how to actually apply eq. (32) for the Wiener Filter, since <span class="math-container">$P'_S$</span> is of course a shorter array than <span class="math-container">$P'_S(k)+P'_N(k)$</span>, 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 <span class="math-container">$\exp{W(k)}$</span>.</p>
<p>When I extrapolate my signal and noise model(and fit only in the appropriate reason) to the full range it looks like this.
<a href="https://i.sstatic.net/4kJTx.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/4kJTx.png" alt="enter image description here" /></a>
The green curve shows <code>P_S+P_N</code>, where I calculated the sum in non-logarithmic space and then transferred back using <code>np.log()</code>. 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 <code>W_log = P_S /(P_S+P_N)</code> the problem is that all my function values are negative, that is why I first use <code>np.exp()</code> 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.</p>
<p><a href="https://i.sstatic.net/cAASP.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/cAASP.png" alt="enter image description here" /></a></p>

CC BY-SA 4.0
2023-05-16T10:34:52.383Z
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
CC BY-SA 4.0
2023-05-16T10:38:34.277Z
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
CC BY-SA 4.0
2023-05-16T11:39:13.837Z
you pointed out the mistake, I mixed up 2 variable names in python. I edited the question with the correct results
trynerror
CC BY-SA 4.0
2023-05-16T12:09:55.110Z
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
<observational-astronomy><galaxy><data-analysis><spectra><python>
How to construct optimal Wiener Filter?
trynerror
# How to construct optimal Wiener Filter?

(Asked by: username_0 on 2023-05-14)

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)}$:
[![enter image description here](https://i.sstatic.net/zHuhh.png)](https://i.sstatic.net/zHuhh.png)

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:

[![enter image description here](https://i.sstatic.net/zAKAr.png)](https://i.sstatic.net/zAKAr.png)

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.
[![enter image description here](https://i.sstatic.net/4kJTx.png)](https://i.sstatic.net/4kJTx.png)
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.

[![enter image description here](https://i.sstatic.net/cAASP.png)](https://i.sstatic.net/cAASP.png)

**Tags:** observational-astronomy, galaxy, data-analysis, spectra, python

> **Comments:**
> **username_1**: 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.  
> **username_1**: 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!)  
> **username_0**: you pointed out the mistake, I mixed up 2 variable names in python. I edited the question with the correct results  
> **username_1**: 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!*  

---

## 2 answers

### Answer by username_1 on 2023-05-14. Score: 2 | Accepted answer

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

> **Comments:**
> **username_0**: 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).  
> **username_0**: 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 ? !  
> **username_1**: @username_0 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?  
> **username_0**: 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  
> **username_1**: @username_0 thanks for the "accept" vote but please consider unaccepting for now until we get this all figured out. It may turn out that you find a bug yourself, in which case it's perfectly okay to post an additional answer to to your own question, and even a few days later accepting it. Goals is always to get things solved :-)  
> **username_0**: I will edit my question in the evening, to show how I resolved the problem.  
> **username_1**: @username_0 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`](https://astronomy.stackexchange.com/posts/53681/revisions)" under the question) and you can copy/paste that to new answer post. *thanks!*  

---

### Answer by username_0 on 2023-05-16. Score: 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.

> **Comments:**
> **username_1**: Yay, congrats on working the problem and finding the solution!  

---

