Extremely low frequency filters
if this is a first-order HPF DC blocking filter, this is the recommendation i have for it in fixed-point, and i think it would be a good idea in floating point, too.

while (n&lt;something)
  {
  y[n]   =  x[n]  -  w[n]*( 2^(1-N) );
  w[n+1] =  w[n]  +  y[n]*( (2^(N-1))*(1-p) );
  n = n+1;
  }


where p is the pole (in the z-plane and inside the unit circle) and $0&lt;1-$p$\ll 1$.  N is a big integer.  with fixed-point you would perform the w[n]*(2^(1-N)) multiplication with arithmetic shift right (ASR) or the (w[n]&gt;&gt;(N-1)) operation.  in float you can just multiply it and that simply changes the exponent in the floating point value.  the constant coefficient (2^(N-1))*(1-p) should be precomputed and N chosen so that this is between $\frac12$ and $1$.

the -3 dB corner frequency $\omega_0 = 2\pi\frac{f_0}{f_\mathrm{s}}$ (as normalized angular frequency) is a function of the pole p and is:

$$\begin{align} 
\cos(\omega_0) &amp;= 1 - \tfrac{(1-p)^2}{2p} \\
 &amp;= 2 - \tfrac12 \left(p + \tfrac1p \right) \\
 &amp;= 2 - \cosh(\ln p)
\end{align}$$

the pole p is

$$\begin{align}
p \ &amp;= \ e^{\operatorname{arccosh}(2-\cos(\omega_0))} \\
\\
 &amp;= \ 2 - \cos(\omega_0) - \sqrt{\big(2 - \cos(\omega_0)\big)^2 - 1} \\
\\
 &amp;= \ 1 - 2\sin\left(\tfrac{\omega_0}{2}\right) \left(\sqrt{1 + \sin^2\left(\tfrac{\omega_0}{2}\right)} - \sin\left(\tfrac{\omega_0}{2}\right) \right) \\
\end{align}$$

so
$$ 0 \ &lt; \ 1-p \ = \ 2\sin\left(\tfrac{\omega_0}{2}\right) \left(\sqrt{1 + \sin^2\left(\tfrac{\omega_0}{2}\right)} - \sin\left(\tfrac{\omega_0}{2}\right)\right) \ \ll \ 1 $$

and the shift or exponent change N is an integer and should be about

$$ N \approx 1-\log_2(1-p) \approx 1 - \log_2(\omega_0) $$

someone can check my math.

this is an improvement by Tim Wescott to my DSP guru trick.  it deals with both the cosine problem and the limitations of fixed point arithmetic.



latest thinking:  i am not sure, but this might work fine for floating-point operation x, y, and w are all single-precision floats (double-precision floats should be no prob).  (1-p) is a coefficient calculated in advance and is very close to zero.

while (n&lt;something)
  {
  y[n]   =  x[n]  -  w[n]*(1-p);
  w[n+1] =  w[n]  +  y[n];
  n = n+1;
  }