Barulho para colorir – Matt Blog

Recentemente, caracterizei um giroscópio de alto desempenho no trabalho por medição Variante de Allan sua taxa de saída. A partir do diagrama de variância de Allan, é possível visualizar e calcular várias fontes de ruído que causam erro, incluindo ruído de quantização, marcha aleatória angular, instabilidade de viés e velocidade aleatória. Cada uma dessas fontes de ruído é causada por um processo físico diferente e pode variar de cor. Qual é a cor do ruído, você pergunta?

A cor do ruído refere-se à forma de sua densidade espectral de potência (PSD). Com o tempo, as pessoas atribuíram cores a certas formas no espectro de potência.

• Branco: Provavelmente a cor de ruído mais famosa, possui um PSD constante para todas as frequências. O fundo “estático” em partes não utilizadas do espectro de rádio está próximo do ruído branco. Em forma de áudio, pode ajudar a bloquear o ruído de fundo ou aliviar os sintomas do zumbido.
• Rosa: Tem um PSD proporcional a $ 1 / f $. Isso significa que no espaço logarítmico ele cai a uma taxa de 10 dB/década. Aparece em muitos componentes eletrônicos como ruído cintilantee os técnicos de áudio costumam usá-lo para caracterizar os alto-falantes.
• Brownovska: Às vezes chamado de ruído “marrom” ou “vermelho”, possui um PSD proporcional a $ 1 / f ^ 2 $. Processos físicos envolvendo caminhada aleatória criam ruído com esta forma PSD.
• Azul: Ruído com PSD aumentando proporcionalmente com $ f $. Os fótons emitidos devido à radiação Cherenkov possuem este PSD. O ruído azul bidimensional é comumente usado em gráficos de computador como um padrão para imagens desfocadas.
• Roxo: Tem um PSD proporcional a $ f ^ 2 $. O ruído acústico na água devido a pequenas mudanças de pressão devido a efeitos térmicos é melhor marcado por essa cor de ruído.

Em última análise, essas cores de ruído geralmente são aproximações do que observamos em certas bandas de frequência, e todo o processo PSD pode conter várias cores de ruído misturadas. Por exemplo, a densidade de ruído da tensão de referência de entrada (que é proporcional à raiz quadrada do PSD) de um amplificador operacional típico é mostrada abaixo. Este ruído pode ser descrito como ruído rosa até cerca de 200-300 Hz e ruído branco até cerca de 50 kHz, mas a densidade final do ruído medido é mais complexa.

Gerando amostras de floresta coloridas

O processo de geração de ruído colorido é relativamente simples:

1. Comece exibindo o sinal de ruído branco no domínio da frequência
2. Formate o sinal no domínio da frequência de acordo com o PSD que você deseja alcançar
3. Tomemos a transformada inversa de Fourier do sinal moldado

Veremos como podemos fazer isso facilmente em Python usando numpy, scipy e matplotlib.

Ruído básico

Para começar, precisamos gerar algum ruído branco no domínio da frequência. Idealmente, esse ruído deve ser normalizado para que tenha um PSD unilateral de 1,0 em todas as frequências, o que simplificará o processo de formatação na próxima etapa.

Existem várias maneiras de fazer isso, e a mais simples é gerar ruído no domínio do tempo e usar a transformada rápida de Fourier (FFT). O PSD de ruído branco discreto amostrado uniformemente é $ 2 \ sigma ^ 2 \ Delta t $, onde $ \ sigma ^ 2 $ é a variância do ruído e $ \ Delta t $ é o período de amostragem. Vamos desligar rapidamente a aula usando esse conhecimento para gerar nosso ruído básico.

import numpy as np


class NoiseGenerator:
    rng = np.random.default_rng()

    def generate(self, dt, n, colour=None):
        """Generates uniformly sampled noise of a particular colour.

        Parameters
        ----------
        dt : float
            Sampling period, in seconds.
        n : int
            Number of samples to generate.
        colour : function, optional
            Colouring function that specifies the PSD at a given frequency. If
            not specified, the noise returned will be white Gaussian noise with
            a PSD of 1.0 across the entire frequency range.
        
        Returns
        -------
        numpy.array
            Array of sampled noise, length `n`.
        """

        if colour:
            raise NotImplementedError("We're not ready yet!")
        
        f, x_f = self._base_noise(dt, n)

        return np.fft.irfft(x_f)

    def _base_noise(self, dt, n):
        """Produces a frequency domain representation of uniformly sampled
        Gaussian white noise.
        """
        
        # Generate Gaussian white noise in the time domain
        sigma = 1 / np.sqrt(2.0 * dt)
        x_t = rng.normal(0.0, sigma, n)

        # Calculate the frequency bins and associated FFT
        f = np.fft.rfftfreq(n, dt)
        x_f = np.fft.rfft(x)

        return f, x_f

Aqui vale a pena verificar se o nosso ruído gerado tem um PSD constante de 1,0. método de Welch ^{é uma técnica comum para estimar sinais PSD no domínio do tempo. Vamos executar algum código para calcular a estimativa PSD e desenhá-la junto com o sinal de ruído.}

import scipy.signal
import matplotlib.pyplot as plt


def plot_psd(x_t, dt, **kwargs):
    """Plots a signal and its PSD given samples and sampling period."""

    t = np.linspace(0.0, dt * len(x_t), len(x_t), endpoint = False)
    f, s_f = scipy.signal.welch(x_t, fs=1/dt, **kwargs)

    plt.subplot(2, 1, 1)
    plt.plot(t, x_t, linewidth=0.5)
    plt.xlabel("Time (s)")
    plt.ylabel("Value (V)")
    plt.subplot(2, 1, 2)
    plt.loglog(f, s_f, linewidth=0.5)
    plt.xlabel("Frequency (Hz)")
    plt.ylabel("PSD ($V^2/Hz$)")
    plt.tight_layout()

# Generate white noise and plot the PSD
dt = 1.0e-3
n = 1000000
noise_generator = NoiseGenerator()
x_t = noise_generator.generate(dt, n)
plot_psd(x_t, dt, nperseg=65536)

A julgar pelo diagrama acima, temos um sinal de ruído branco com PSD próximo a 1,0 na maior parte da banda de frequência.

Quem assistiu a implementação NoiseGenerator Eu observaria cuidadosamente a ineficiência de gerar ruído aleatório no domínio do tempo, depois transformar no domínio da frequência e depois reconverter para o domínio do tempo. É possível gerar ruído branco diretamente no domínio da frequência^{. Aqui está uma implementação alternativa _base_noise() ele faz exatamente isso.}

    def _base_noise(self, dt, n):
        # Calculate random frequency components
        f = np.fft.rfftfreq(n, dt)
        x_f = self.rng.normal(0,
            0.5, len(f)) + 1j * self.rng.normal(0, 0.5, len(f))
        x_f *= np.sqrt(n / dt)

        # Ensure our 0 Hz and Nyquist components are purely real
        x_f[0] = np.abs(x_f[0])
        if len(f) % 2 == 0:
            x_f[-1] = np.abs(x_f[-1])
        
        return f, x_f

Projeto de espectro

Neste ponto, tudo o que resta é colorir o ruído moldando seu espectro. Atualizar generate() manipular a função de coloração é uma maneira genérica apropriada de fazê-lo. Adicionaremos algumas funções estáticas auxiliares à classe para as cores usuais, mas, claro, você pode fazer algo completamente personalizado. Há realmente apenas duas coisas para mencionar aqui:

1. O PSD é diretamente proporcional ao quadrado da FFT, então precisamos tirar a raiz quadrada de nossas funções de coloração.
2. Várias cores de floresta têm $ f $ no denominador da fração, portanto, precisamos abordar especificamente o caso em que $ f = 0 $.

    def generate(self, dt, n, colour=None):
        f, x_f = self._base_noise(dt, n)

        if colour:
            x_f *= np.sqrt(colour(f))

        return np.fft.irfft(x_f)
    
    @staticmethod
    def pink(scale=1.0):
        """Creates a pink noise colouring function.

        Parameters
        ----------
        scale : float, optional
            Multiplier to adjust the scale of the pink noise.
        
        Returns
        -------
        function
            Pink noise colouring function that can be used with `generate()`.
            The function returned will be :math:`y(f) = s / f`, where :math:`s`
            is `scale`. At f = 0, the function will simply return 0.0.
        """

        return lambda f: scale / np.where(f == 0.0, np.inf, f)

Deixei de fora a maioria das cores no trecho de código acima, mas o todo NoiseGenerator classe pode ser encontrada em essa essência. É de especial interesse função de colorir por peçasque permite a criação de algum ruído funky que não pode ser classificado por nenhuma das cores padrão.

ng = NoiseGenerator()
noise_settings = {
    "Piecewise": {
        "dt": 1e-4,
        "n": 1000000,
        "colour": ng.piecewise_logarithmic(
            [90.0, 100.0, 110.0, 450.0, 500.0, 550.0],
            [0.1, 10.0, 0.01, 0.01, 2.0, 0.001]
        )
    },
    "Brownian": {"dt": 1e-3, "n": 100000, "colour": ng.brownian()},
    "Azure": {"dt": 1e-3, "n": 100000, "colour": ng.azure(1.0e-6)}
}
for name, settings in noise_settings.items():
    x_t = ng.generate(**settings)
    plot_psd(x_t, settings["dt"], nperseg=65536)
plt.legend(noise_settings.keys())

Perguntas feitas

Embora a classe de geração de ruído de cor atendesse perfeitamente às minhas necessidades, algumas dúvidas surgiram enquanto eu a usava. Por enquanto, deixei-os sem resposta, mas se alguém tiver respostas, avise-os.

Distribuição de ruído

A função de densidade de probabilidade usada na geração anterior de números aleatórios é gaussiana. É importante notar que a distribuição de ruído é um conceito separado inteiramente em relação à sua densidade espectral de potência. Para nossa geração de ruído de tempo base, poderíamos usar uma distribuição uniforme, ou uma distribuição exponencial, ou uma distribuição gama, e tudo funcionaria.

Curiosamente, após a transformada de Fourier, todas essas distribuições produzem algo que se parece aproximadamente com uma distribuição gaussiana 2D em um plano complexo. Além disso, depois de modelar o PSD usando a função de coloração e tomando a transformada inversa de Fourier, sempre parece que acabamos com ruído no domínio do tempo que tende à distribuição gaussiana.

1. É possível gerar ruído com uma distribuição específica no domínio da frequência?
2. Por que o ruído de formatação no domínio da frequência tende a deslocar a distribuição no domínio do tempo em direção à distribuição gaussiana? Isso tem alguma coisa a ver teorema da fronteira central?
3. É possível gerar ruído com uma determinada distribuição e uma determinada cor que não seja branca?

Integração de ruído

Usei ruído gerado sinteticamente para simular um giroscópio, então integraria um sinal para calcular o erro angular do sistema. Dessa forma, pude testar algoritmos de estimativa de viés e compensação sem a necessidade de coletar longos conjuntos de dados de um giroscópio real. Ainda assim, notei um comportamento estranho no sinal integrado. Para todo ruído colorido que tivesse componente de frequência de 0 Hz 0,0, a integral resultante teria um valor inicial e final de 0,0. Obviamente, isso significa que, se eu gerar 1 milhão de amostras de ruído, meu erro angular por milionésima amostra será sempre 0,0.

O que está acontecendo aqui? Minha teoria atual é que isso se deve à definição de uma transformada discreta de Fourier:

$$ X_k = \ soma_ {n = 0} ^ {N-1} x_n \ cdot e ^ {- \ frac {i 2 \ pi} {N} kn} $$

Se a transformada de Fourier for 0,0 para $ k = $ 0, o acima simplifica para:

$$ 0,0 = \ soma_ {n = 0} ^ {N-1} x_n $$

Claro, isso é exatamente o que vemos incomum – a soma de nossas amostras de ruído é 0,0. Isso levanta a questão de se definir um componente de 0 Hz em nossas funções de modelagem de cores para 0,0 é razoável. Talvez devesse ser algum outro valor?

Observações finais

Este pequeno mergulho nas técnicas de geração de ruído resultou em um NoiseGenerator uma classe que tem sido útil para simulações de giroscópio, bem como outras aplicações. Com apenas 44 linhas de código funcional, é uma prova do poder das bibliotecas de computadores científicos no ecossistema Python. Embora a geração e filtragem de ruído branco pareçam ser uma técnica aceita para gerar ruído colorido, levantou algumas questões sobre a estrutura do ruído gerado que espero sejam respondidas no futuro.