How to Create a Filter that Matches a Desired Magnitude Response

Question

I'm very, very new to this stuff, so please forgive me if the answer is obvious. I haven't been able to find any clear answers to this issue, even after looking through source code for various filter implementations and scouring the net for similar questions.

I have a function $f(x)=y$, such that:

• $x$ is a frequency in Hz
• $y$ is an output dB it should, ideally, be quieted down by.

I would like to apply this function to a digitally-stored sound - represented of course by samples in the time domain, as is standard.

Ideally, I would like this to give the same output as if one manually took each frequency out of a sound, amplified them with the function's output, then added them back together. I know this method itself is beyond impractical, so I'm sure my answer lies somewhere in the realm of IIR or FIR filtering algorithms.

As an aside, I would prefer to avoid approximatory approaches where possible; I'd like to get as close as I can to applying the function exactly, rather than approximating its curve. This is for research purposes - I don't plan to apply this filter in real time, and my programming language of choice (Lua) isn't really suited to that anyways - so the increased speed and memory usage this might incur is not of issue to me.

Where should I begin? Is this possible to do?

Answer 1

What you seem to want is a multiplicative modification of the spectrum of the input signal in the following form:

$$Y(\omega)=H(\omega)X(\omega)\tag{1}$$

where $X(\omega)$ is the spectrum of the input signal, $Y(\omega)$ is the spectrum of the output signal, and $H(\omega)$ is the function that modifies the input spectrum.

Equation $(1)$ is exactly what a (linear and time-invariant) filter does. Examples of such filters in the discrete domain are FIR and IIR filters. In general, the function $H(\omega)$ - called the filter's frequency response - is complex-valued. That means that not only the magnitude of the input spectrum is changed but also its phase. If I understand your question correctly, you just care about the magnitude.

For clarification, the frequency response's influence on a sinusoidal input signal $x[n]=\sin(\omega_0n)$ is described by the following equation:

$$y[n]=\big|H(\omega_0)\big|\sin\big(\omega_0n+\phi(\omega_0)\big)\tag{2}$$

where $\phi(\omega)$ is the phase of the complex-valued frequency response $H(\omega)$:

$$H(\omega)=\big|H(\omega)\big|e^{j\phi(\omega)}$$

Eq. $(2)$ shows that the magnitude of a sinusoid at frequency $\omega_0$ is modified by the magnitude of the frequency response at that frequency. Furthermore, its phase is modified by the phase of the frequency response at that frequency.

The catch is that you can't just prescribe some $H(\omega)$ and expect a realizable filter to exactly implement the given response. You will always have to live with an approximation of the desired response. However, if computational cost and memory requirements are of no concern you can get very close to the prescribed response.

I would recommend to use an FIR filter to approximate the given desired response. Reasons why one might want to avoid FIR filters are their large delay (at least for linear phase filters) and their large computational complexity, both when compared to equivalent IIR filters. However, I think that for your application these drawbacks are not very relevant, and the advantages of FIR filters outweigh the disadvantages. The advantages are ease of design and inherent stability.

The frequency response of a causal FIR filter is

$$H(\omega)=\sum_{n=0}^{N-1}h[n]e^{-jn\omega}\tag{3}$$

where $h[n]$ are the filter coefficients ("taps"), and $N$ is the filter length (= number of taps). Obviously, the larger $N$, i.e., the more taps, the better the approximation of a given response will be.

Given the filter coefficients $h[n]$, the output $y[n]$ is computed by the discrete convolution of the filter coefficients with the input $x[n]$:

$$y[n]=\sum_{k=0}^{N-1}h[k]x[n-k]$$

Designing a filter means finding the "best" filter coefficients $h[n]$ for a given specification. There are countless filter design methods and it's not always straightforward to find the one that best fits a given problem. In your case I would probably start experimenting with the frequency sampling method. In this case, the desired response is sampled on a dense equidistant frequency grid. This discretized response is then transformed to the time domain by applying an inverse discrete Fourier transform (IDFT), which is usually implemented by the (inverse) Fast Fourier Transform (IFFT). Finally, the resulting impulse response is usually windowed to smooth out the resulting response and to reduce the impulse response to the desired filter length $N$.

For illustration purposes I used the frequency sampling method to design a filter approximating a rather silly fantasy specification. The desired magnitude response has a staircase shape and was sampled at $1000$ equidistant points. The desired filter length was chosen to be $N=201$. The figure below shows the result. Clearly, an actual filter cannot perfectly follow the discontinuities of the desired response. However, increasing the filter length will make the transitions narrower.

Here is a link to a simple Octave/Matlab implementation of the frequency sampling method used in the design example shown above. The filter was designed with the following commands:

mag = [ones(1,200),.75*ones(1,200),.5*ones(1,200),.25*ones(1,200),zeros(1,200)];
h = fsamp( mag, 201 );

Answer 2

If I understood the question correctly, here's one possible method to manipulate audio signal in sample level.Method is based on one pole IIR filter which actually only changes the value of each sample as wanted. If one wants to place target values for output samples then some additional basic math is needed.

% Octave packages -------------------------------
pkg load signal

clear all;
clf;

fs = 1000;
fc = 50;

% random signal
N = fs;
x = 0.25*randn(N,1); %random signal 
%x = [ones(1,200),.75*ones(1,200),.5*ones(1,200),.25*ones(1,200),zeros(1,200)];
y = zeros(N,1);

% list of frequencies and gains for manipulation
fcoctaves = [fc, fc*2, fc*3, fc*4, fc*5]; %Hz
gainsoct = [-3, -6, -12, -24, -32]; % dB

% output global gain change (dB)
g = 0;

% global gain change loop
for n = 1:1:fs
  y(n) = x(n) * dB2lin(g);  % DF1, b(1) = dB2lin(g)
endfor

% fixed frequency gain changes
for n = 1:length(fcoctaves)
    y(fcoctaves(n)) = x(fcoctaves(n)) * dB2lin(gainsoct(n)); % DF1
endfor

audiowrite('output.wav', y, fs);

figure(1);
plot(x, "linewidth", 1, '-x-')
hold on;
plot(y,"linewidth", 1 ,'--o')
legend('Original response', 'Modified response',  'location', 'southwest');
grid on;

dB2lin function:

function linampl = dB2lin(dB)
  linampl = 10^(dB/20);
endfunction

Example source code handles only one second of audio data but it's not difficult to handle usage with longer audio data as one gets all needed info from source audio file (fs, length). Example code results:

I tried also example data Matt gave in his answer with additional manipulation of few samples as in my example. Result:

As I don't have much experience in DSP I can't say if this works as what OP asked nor if there are bottlenecks turning this method unusable.