How to create a filter that is a single period of a time-warped raised cosine?

Question

I am attempting to implement a filter described in the paper Functional Count-Comparison Model for Binaural Decoding by Pulkki et al 2009:

Here is Figure 9:

However, I am quite certain that the equation for n is incorrect since the cosine is squared and there is a fourth root of n/fs. This equation does not produce a single period for every omega value (in my case between 100 Hz and 12 kHz).

I am also unsure why fs is in the equation for f at all, since it should be handled by the length of n. Am I correct on either or both of these points?

Another issue I am having is that the BF of the highest frequency band is 12 kHz, while fs = 20 kHz, which I suspect will cause some aliasing. Should I remove all frequency bands above 10 kHz?

The final issue I have is that given the equation, there does not seem to be a high enough fs to represent the filter for the highest frequencies as the length of n quickly goes to 0 well before 10 kHz.

Here is my attempt at implementing it in python:

import numpy as np
import matplotlib.pyplot as plt


# Define the function
def raised_cosine_filter(omega, fs=20000):
    # Calculate the period based on omega
    period = (2 * np.pi) / omega

    # Calculate the number of samples
    num_samples = int(np.ceil(period * fs))

    # Generate x values for one period
    n = np.linspace(0, period, num_samples)

    return ((np.cos(omega * (n / fs) ** 0.25 - np.pi) + 1) ** 2)


# Example usage:
f = 100  # Frequency in Hz
omega = 2 * np.pi * f  # Angular frequency in rad/s
y = raised_cosine_filter(omega)

# Plot the raised cosine pulse
plt.plot(y)
plt.xlabel('Sample number')
plt.ylabel('y')
plt.title(f'One Period of the Raised Cosine Filter (f = {f} Hz)')
plt.grid()
plt.show()

And here is the resulting plot, which looks incorrect as you can see:

Any help would be appreciated in fixing either my implementation or clarifying whether the paper has a mistake. Thank you.

Answer 1

According to the erratum

Equation 5 is an error. Here is the actual equation:

And the code after the correction:

import numpy as np
import matplotlib.pyplot as plt


# Define the function
def raised_cosine_filter(fbf, fs=20000):

    # Calculate the number of samples
    num_samples = int(np.ceil(fs / fbf))

    # Generate x values for one period
    n = np.arange(num_samples)

    return 0.25 * (np.cos(2 * np.pi * (n * fbf / fs) ** 0.25 - np.pi) + 1) ** 2


# Example usage:
fbf = 750  # Frequency in Hz
y = raised_cosine_filter(fbf)

# Plot the raised cosine pulse
plt.plot(y)
plt.xlabel('Sample number')
plt.ylabel('y')
plt.title(f'One Period of the Raised Cosine Filter (f = {fbf} Hz)')
plt.grid()
plt.show()

Here is what the plot looks like:

I'm still concerned about the sampling rate, but not sure how to address that without making the filter inconsistent with the data I'm using.

Answer 2

Not a definite answer, but a tentative hint that's too long to put in a comment

---

I could be completely off mark but I believe there's either a typo in the paper or no typo but rather something off with your implementation. Regardless, if you take $f_s$ out of the equation for $f$ like you hinted at, and use *normalized* angular frequency, i.e. $\omega = 2\pi f/f_s$, you get something that looks very close (up to normalization) to what you're looking for.

I used $750 \,\texttt{Hz}$ to match what the authors do in the paper linked:

import numpy as np
import matplotlib.pyplot as plt


# Define the function
def raised_cosine_filter(omega, fs):
    # Calculate the period based on omega
    period = (2 * np.pi * fs) / omega

    n = np.linspace(0, period, 360)

    return (np.cos(omega * n ** 0.25 - np.pi) + 1) ** 2


# Example usage:
f = 750  # Frequency in Hz
fs = 20000 # sampling frequency
omega = 2 * np.pi * f / fs  # normalized angular frequency 
y = raised_cosine_filter(omega, fs)

# Plot the raised cosine pulse
plt.plot(y)
plt.xlabel('Sample number')
plt.ylabel('y')
plt.title(f'One Period of the Raised Cosine Filter (f = {f} Hz)')
plt.grid()
plt.show()