That are just the same questions I’ve been pondering over. Brown noise should describe a “Brownian” motion. I’ve not looked at the concrete physical conditions that make up this kind of noise.
Perhaps it would be best to start the input signal with a zero in first place. But as you’ve said it is theoretically possible to throw 1000 times a 6. Thus we would stay in the positive region for a long time. Folding back seems a neat solution, especially if you want to preserve the initial value.
The other thing is to place the extreme peaks equally away from 0 and to normalise the signal - which is of course nothing but removing the DC offset. However,the signal could still reside in the positive/negative region for a long time, as you’ve remarked rightly.
It is very difficult to translate the laws of infinite signals to discreet signals of finite length.
You’ve raised an interesting question, namely: How to construct a LF control signal from our noise functions?
If I want a signal from 1 to 20 Hz, this would mean that the amplitude at 1 Hz is 0 dB and at 16 Hz already -24 dB. However well constructed our noise may be, the application of the highpass (for DC) and perhaps a lowpass at 20 Hz will utterly destroy the nicely constructed frequency response, at least I’ve got this suspicion.
what would be the best to do for the above mentioned signal? Should we lower the cut-off frequency or should we take the signal from a 20 times oversampled noise?
Maybe it is best to test our filter firstly with sine waves for the corresponding frequencies (1 2 4 8 10 16…) and to see how the amplitudes evolve. You can of course simply take a look at the spectrum plot (although it is difficult to judge within this regions).
We can remove DC without filtering:
(defun remove-dc (sig)
(diff sig (calc-dc sig)))
;;; calculate dc offset
(defun calc-dc (sig)
(sref (sound (mult (integrate sig)
(/ (get-duration 1))))
(/ (1- len) *sound-srate*)))
(remove-dc s)
To limit the upper frequency range we can make use of the Nyquist–Shannon theorem and the “resample” function
If we resample noise to a low sample rate, we know that the maximum frequency will be limited to half the sampling rate.
(defun lowpass (sig hz)
(resample
(force-srate (* hz 2) sig)
(snd-srate sig)))
;; lowpass filter at 20 Hz
(lowpass s 20)
Putting this all together in an example, for brown noise, band-limited from DC to 100 Hz:
;;; DC removal
(defun remove-dc (sig)
(diff sig (calc-dc sig)))
;;; calculate dc offset
(defun calc-dc (sig)
(sref (sound (mult (integrate sig)
(/ (get-duration 1))))
(/ (1- len) *sound-srate*)))
;;; extreme slope low-pass filter
(defun lowpass (sig hz)
(resample
(force-srate (* hz 2) sig)
(snd-srate sig)))
;;; brown noise generator
(defun brown ()
(mult (/ 200 (sqrt (get-duration 1))) ; normalization fudge factor
(integrate (noise))))
(remove-dc (lowpass (brown) 100))
Hi Steve
I am afraid your code won’t work like that. There are several problems we’re running into.
The first one - which is quite subtle and easily correctable - is the calculation of the DC-offset.
As I’ve mentioned in an earlier thread, the integration of a sound does not return the last sum. This can blur the result evidently. Especially in the case of our brown noise, where the difference of the last two samples can be up to 2.
To get the last summation, at least one 0 must be added to the signal.
However, the same result (the mean or DC-offset) can be got with a little less of code by simply averaging the signal.
;; a sample sound to illustrate the problem
;; Can be omitted to see the results
;; for the actual selection
(setf s (seq (snd-const -0.1 0 *sound-srate*
(/ (1- len) *sound-srate*))
(snd-from-array 0 *sound-srate* #(1))))
;; DC as simple average
(setf ln1 (truncate len))
(print (snd-fetch (snd-avg s ln1 ln1 op-average)))
;; DC from running sum without correction
(print (sref (sound (mult (integrate s)
(/ (get-duration 1))))
(/ (1- len) *sound-srate*)))
;; DC from running sum with added silence for correction
(print (sref (sound (mult (integrate (seq s (s-rest 1)))
(/ (get-duration 1))))
(/ len *sound-srate*)))
One can easily see how the second value deviates from the first and last one.
Normally, when the sound is quite long, this is no problem. However, it is significant when we want to return the average RMS value of a sound, where 4410 samples (10 Hz) could be hidden within the very last sample value of the RMS signal.
It is possible that your method is faster than mine and less resource-killing. Snd-avg could easily performed in a recursive manner to free memory between stages.
The DC-problem could also be approached with a simple DC-blocker.
(defun dc-blocker (s-in 3db-freq)
(let* ((pole (- 1 (/ (* pi 2.0 3db-freq) *sound-srate*))))
(biquad-m s-in 1 -1 0 1 pole 0)))
(dc-blocker s 40.0)
This has the advantage that DC level-changes are followed up.
However, a “absolute” correction maybe still necessary (the filter needs some time to establish; you can try the code before the print commands above to test it) and it is for sure not ideal when we want to use the low frequencies as a control signal (although the attenuation is not too strong). I fancy, it is rather a method for real time manipulations.
The second problem is much more severe. Your low pass won’t work as such. Due to the Nyquist/Shannon theorem, the higher frequencies will indeed be blocked and …
… fold back in the audible audio band.
;;lowpass10 and test
(defmacro sweep (dur)
`(abs-env (mult (pwl 0.0 0.8 ,dur 0.05 ,dur)
(hzosc (sum 440 (mult 1760 (ramp ,dur)))))))
(defun lowpass (sig hz)
(resample
(force-srate (* hz 2) sig)
(snd-srate sig)))
(defun lowpass10 (x hz)
(lowpass2 (lowpass2 (lowpass2 (lowpass2 (lowpass2 x
hz 3.1969)
hz 1.1013)
hz 0.7071)
hz 0.5612)
hz 0.5062))
(seq (sweep 3)
(lowpass10 (sweep 3) 880)
(lowpass (sweep 3) 880))
Three sweeps are produced, The first is not damped, the second with a 10th-order butterworth and the last by resampling. The frequencies are such that they can easily be heard.
It is a good way too illustrate aliasing.
I would have expected that the resample function (which is high-quality-interpolation) could deal too some degree with the aliasing (of course only if it is used for downsampling too).
But the function seems to be limited, it only returns some crackle when used for down sampling (a bug?). At higher frequencies it works however.
Fortunately, we can produce the brown noise from the very start at a lower samplerate and then resample it.
If the brown noise should be returned to Audacity, it is even simpler because all sounds are automatically seen as if they were at the track sample rate.
For example: if we produce a sine wave at 440 Hz and a sample rate of 192000 Hz, Audacity will but put the incoming samples into the track (which has for simplicity 48000 Hz). Thus the sound will be 4 times longer and the new frequency is therefore 110 Hz.
Thus, to produce brown noise that starts at 1 Hz (instead of 20 Hz), we have to oversample it by a factor of 20 (and a duration of 1/20) and pass it back to Audacity.
Admittedly, the problem of the higher frequencies is still here, if we want a smoother curve, we still have too employ a lowpass filter. The question is if the brown noise is not already too much attenuated with the resampling method.
The start of the 6 dB roll-off seems to be widely accepted especially in by synthesizer manufacturers and users.
It is obvious that this limit of 20 Hz is moved down by our method, such that the sound will already be attenuated by 30 dB at 32 Hz. It is a question of trial and error, if a further damping is needed. One has too bare in mind that Lowpass filters that are implemented as IIR (such as biquads) do not work very well under 1 kHz due to quantisation noise (accumulated in the feedback loop) and numerical round-off errors. We simply replace sophisticatedly created noise with system created one.
A FIR filter would work much better (with the advantage of phase linearity above all) but is very slow (with increasing steepness).
You’re quite right, after the “quick” method there is still a liitle residual DC offset, but only about -104 dB, which for most purposes is probably not very significant. Yes, if you need it more accurate you can pad the sound with an extra sample. Alternatively, to keep the error extremely small, but still reasonably efficient you could use a hybrid:
;;; calculate dc offset
(defun calc-dc (sig)
(if (< (snd-length sig 1000000) 1000000)
(snd-fetch
(snd-avg s (truncate len)(truncate len) op-average))
(sref (sound (mult (integrate sig)
(/ (get-duration 1))))
(/ (1- len) *sound-srate*))))
Snd-avg can only be used to a certain degree without several stages, i.e. when len is not too big. Therefore, you’ll better use the integreation with the 0 sample tail.
Or you use the following neat replacement for the integrate function, where the multiplication with sound-srate can be dropped and the resulting sound’s length is always equal to the input sound. In other words, all summations are made.
(defun integrator (s )
(biquad s 1 0 0 1 1 0))
(setf s (hzosc 0.9))
(print (/ (sref (sound (integrator s))
(/ (1- len) *sound-srate*)) len))
(print (sref (sound (mult (integrate
(seq s (snd-from-array 0 *sound-srate* #(0))))
(/ (get-duration 1)))) (/ len *sound-srate*)))
(setf ln1 (truncate len))
(print (snd-fetch (snd-avg s ln1 ln1 op-average)))
The biquad solution is neat, but no solution will be as quick as just looping through the samples in C, but whichever method is used the Nyquist primitives loop through all of the samples. The best that we can do is use an “efficient for Nyquist” method.
I think that our examples are quite fast for Nyquist code that verges in such a manner to DSP implementations.
Most functions are very fast as long as we do not try to gather information.
For example, sref is very slow because it evalueates the whole sound to get the last value. There is unfortunately no alternative. I think I’ve tried every imagineable method.
- resampling
- indexing with snd-compose
- extracting
- concatinating with a time-shifted sound
- destructively removing of blocks with snd-fetch
- …
Each method works and each method is equally slow.
It’s a pity that sounds can only be accessed sequentially.
A relative or blockwise access method would be a great improvement.