3
$\begingroup$

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?

$\endgroup$
4
  • 1
    $\begingroup$ To confirm, you seek a digital filter that would best match a target magnitude response, and you have no concern with the phase response. There will be no avoiding approximations since an exact match to any arbitrary magnitude function is not possible (in particular step changes in magnitude), but depending on how much filter complexity and delay is tolerable, the match can get to within an allowable error under a specified criteria (such as minimum least squares, for example). Luckily one of our leading contributors wrote his PhD on this, so I believe you'll see a very good answer shortly. $\endgroup$ Mar 10, 2023 at 3:16
  • $\begingroup$ In the meantime check out: mattsdsp.blogspot.com/p/phd-thesis.html which covers simultaneous magnitude and phase constraints but MattL can likely direct you to best approaches for magnitude only concerns. $\endgroup$ Mar 10, 2023 at 12:44
  • $\begingroup$ If you're constraining yourself to Lua -- switch to Python. Lua is pretty close to a subset of Python, and Python + numpy + scipy has become a serious language of choice for a number of signal processing applications. If you don't mind folks who know both languages say "goodness, this Python code looks like Lua" you'll be able to pick up Python pretty quickly, and you won't be writing the signal processing math by hand. If you're in an environment that constrains you to Lua -- you have my sympathies. Prototype your algorithms in Python, then translate to Lua. $\endgroup$
    – TimWescott
    Mar 10, 2023 at 22:40
  • $\begingroup$ If this is for audio -- i.e., if you're trying to equalize sound files for listeners, then both the phase and the amplitude do matter (search on "psychoacoustics" -- filters with different phase responses will sound different, even if they have the same amplitude response). For that matter, phase also matters in some physical systems analysis and in closed-loop control. So you may want to edit your question to make it clear whether you only care about amplitude, or whether phase might matter. $\endgroup$
    – TimWescott
    Mar 10, 2023 at 22:46

2 Answers 2

5
$\begingroup$

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.

enter image description here

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 );
$\endgroup$
0
0
$\begingroup$

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:

enter image description here

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

enter image description here

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.