Barulho para colorir – Matt Blog

By | Junho 16, 2022

Gerando ruído colorido para simular processos físicos


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?

Gráficos escuros com 5 séries de dados para roxo, azul, branco, rosa e marrom, cada um colorido em sua própria cor.
Densidade de potência espectral de algumas das cores florestais mais comuns

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.

Um diagrama de ruído típico de uma folha de dados opamp mostrando duas curvas que diminuem linearmente no espaço log-log até o ponto em que se tornam horizontais.  As duas curvas representam a densidade de ruído de tensão para diferentes tensões de alimentação.
Tensão de referência de densidade de ruído para ADA4625-1 opamp

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)

Uma imagem com dois sublinhados contendo um traço azul cada.  A parte superior exibe um sinal de ruído baseado em um tempo de cerca de 1000 segundos e com picos de até 100 V ou mais.  A parte inferior é um diagrama logarítmico que mostra uma estimativa da densidade de densidade de potência espectral que é relativamente próxima de 1,0 entre frequências de 0,001 a 500 Hz.
Sinal de ruído e sua densidade espectral de potência

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())

Diagrama duplo mostrando 3 sinais de ruído no domínio do tempo e sua densidade espectral de potência.  Ruído de peça com dois picos diferentes em PSD, ruído browniano com diminuição acentuada de 20 dB/década e ruído azure aumentando em PSD com valores muito baixos.
Alguns exemplos de ruído colorido, incluindo ruído com densidade de potência espectral ajustada por peça

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?

Seis parcelas privadas.  No canto superior esquerdo, um diagrama com dois traços - um mostrando ruído distribuído gama (que tem todas as amostras maiores que 0) e outro mostrando ruído rosa uniformemente distribuído em cada lado do eixo x.  No canto inferior esquerdo, os PSDs padrão parecem ruído branco e rosa.  No canto superior direito estão dois histogramas no domínio do tempo - um mostrando a distribuição gama e outro mostrando o que parece ser uma distribuição gaussiana.  No canto inferior direito estão histogramas de domínio de frequência 2D que parecem aproximadamente gaussianos.
À esquerda, gráficos de clima e PSD de ruído branco básico com distribuição gama e ruído rosa formados a partir desse ruído básico. À direita estão os histogramas dos domínios do tempo e da frequência do mesmo. Histogramas de domínio de frequência estão localizados em um plano complexo.

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.

Um diagrama de seis sinais semelhantes a caminhadas aleatórias que começam com um valor de 0,0 e terminam com um valor de 0,0, mas não são correlacionados.
Integral (soma) de seis conjuntos de 1 milhão de pontos de ruído rosa gerado com um PSD de 0,0 a 0 Hz

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.

Deixe uma resposta

O seu endereço de email não será publicado.