Funções de janela DIY

No processamento de sinal digital, as funções de janela são amplamente utilizadas para limitar o sinal no tempo, e seus nomes são bem conhecidos por todos que já encontraram a transformada discreta de Fourier de uma forma ou de outra: Hann, Hamming, Blackman, Harris e outros. Mas são suficientes, é possível inventar algo novo e isso tem algum sentido?



Neste artigo, veremos a saída de uma função de janela com novas propriedades usando o Wolfram Mathematica. Também presume-se que o leitor tenha uma compreensão geral do processamento de sinais digitais no contexto da questão em discussão e esteja pelo menos familiarizado com o artigo da Wikipedia .







Introdução



Historicamente, as primeiras funções de janela surgiram no processo de tentativas de melhorar suas propriedades espectrais pela redução da amplitude dos lobos laterais - já que quando o sinal é multiplicado pela função de janela, seus espectros são convolutos, o que introduz algumas limitações para análise. Naquela época, década de 50 do século 20, o poder de computação não permitia manipular de forma fácil e natural os cálculos simbólicos, e a seleção dos parâmetros ótimos era bastante difícil. Essa é uma das razões pelas quais as funções recebem nomes diferentes - elas, ou melhor, suas propriedades espectrais, já há muito tempo são estudadas por diferentes pesquisadores. Um efeito colateral disso é que o conjunto existente de funções de janela nomeadas é percebido como uma espécie de "conjunto padrão" fora do qual nada pode ser encontrado;ao mesmo tempo, os nomes dessas janelas não trazem nenhuma informação sobre as propriedades e sugerem seu estudo e memorização separados.



Classificação



Todas as funções de janela são obtidas multiplicando alguma função por uma retangular, o que é garantido para garantir sua restrição de tempo. Portanto, em primeiro lugar, eles podem ser classificados de acordo com as funções que utilizam:



  1. a soma dos cossenos. A classe mais extensa devido ao fato de que seu espectro é fácil de calcular e é a soma das funções sinc ponderadas. Isso inclui Hann, Hamming, Blackman, Harris ...
  2. polinômio contínuo por partes. Eles são obtidos como resultado da convolução de funções simples - por exemplo, retangular e triangular. Ao mesmo tempo, seus espectros são multiplicados e sua descoberta também não é particularmente difícil. Estes incluem Dirichlet, Triangular, Parzen, Welch ...
  3. todos os outros - usando exponenciais, Gaussianos, sinc e outros, a escolha de funções específicas em que é mais ideológica por natureza do que propriedades espectrais especificamente.


Eles também podem ser divididos por propriedades nas bordas:



  1. sem lacunas nas bordas - valores iguais a zero. As pausas são de Dirichlet, Hamming, Blackman-Harris. Não tem - Triangular, Hann, Nuttal
  2. ausência de descontinuidades do 1º e outros derivados


Mais duas funções de janela se destacam separadamente:



  1. Kaiser, que permite definir a largura da pétala principal;
  2. Dolph-Chebyshev, todos os lóbulos laterais iguais a uma dada amplitude.


Ambos têm descontinuidades nas bordas dos valores e derivados, e seus cálculos envolvem algumas dificuldades - a função Kaiser é calculada por meio de uma função especial (Bessel) e a função Dolph-Chebyshev é determinada no domínio da frequência e é calculada por meio da transformada discreta inversa de Fourier. Encontrar suas antiderivadas analíticas também é particularmente difícil.



Você pode examinar as funções da janela para intervalos com o seguinte código simples:
wins = {DirichletWindow, HammingWindow, BlackmanWindow, 
   BartlettHannWindow, BartlettWindow, BlackmanHarrisWindow, 
   BlackmanNuttallWindow, BohmanWindow, ExactBlackmanWindow, 
   FlatTopWindow, KaiserBesselWindow, LanczosWindow, NuttallWindow};
Table[{f, 
  Table[Limit[D[f[x], {x, k}], x -> 1/2, 
    Direction -> "FromBelow"], {k, 0, 4}]}, {f, wins}]</code>
<img src="https://habrastorage.org/webt/rm/wq/8f/rmwq8fjkjxjcox-czh8zevkbd0m.png" />

<code>Plot[Evaluate[Through[wins[x]]], {x, -1, 1}, PlotRange -> {-0.25, 1}, 
 GridLines -> {{-0.5, -0.25, 0.25, 0.5, 1}}, AspectRatio -> 5/8, 
 PlotLegends -> "Expressions"]








Engenharia reversa



Vejamos as definições das funções Blackman e Nuttel:



BlackmanWindow[x] // FunctionExpand




{150(25cos(2πx)+4cos(4πx)+21)12x120|x|>12



NuttallWindow[x] // FunctionExpand




{121849cos(2πx)+36058cos(4πx)+3151cos(6πx)+8894225000012x120|x|>12



Eles são somas de ondas de 3 e 4 cossenos com frequências pares (começando do zero) e alguns coeficientes. De onde vêm esses coeficientes? É improvável que eles tenham sido selecionados à mão, pelo menos cada um separadamente - afinal, existem condições de contorno que as janelas devem cumprir independentemente de suas propriedades espectrais - pelo menos iguais à unidade no centro.



Vamos tentar "inventar" a função Blackman por conta própria. Para fazer isso, definimos uma função com coeficientes ainda desconhecidos



f = Function[x, a + b Cos[2 Pi x] + c Cos[4 Pi x]];


e compor um sistema de equações definindo as condições de contorno - igualdade à unidade no centro, zero nas bordas e igualdade a zero das derivadas nas bordas:



ss = Solve[{
   f[0] == 1,
   f[1/2] == 0,
   f'[1/2] == 0
   }, {a, b, c}]








{{b12,c12a}}



A solução foi encontrada, mas não uma, mas muitas - dependendo do coeficiente a , sobre o qual Wolfram educadamente nos alertou. Agora, a partir da solução encontrada, definimos uma nova função dependendo do coeficiente desconhecido: É fácil ver que para a = 0,42 obtemos a função de Blackman. Mas por que exatamente 0,42? Para responder a essa pergunta, precisamos traçar seu espectro. A transformação analítica de Fourier não é a tarefa mais fácil, mas o Wolfram também pode lidar com isso, nos poupando da tarefa.



hx = Function[{x, a},

Piecewise[{{f[x] /. ss[[1]], -1/2 < x < 1/2}}, 0] // Evaluate]
















hw = Function[{w, a}, 
  FourierTransform[hx[x, a], x, w] // #/Limit[#, w -> 0] & // 
    Simplify // Evaluate]






Nota
, FourierTransform , — , , — . w. — w, , .







Desta forma, o gráfico da função ainda não é muito informativo, pois é mais conveniente analisar o espectro em escala logarítmica. Adicionando a conversão apropriada para decibéis , obtemos o mesmo gráfico usual do espectro de funções de janela:20log10(|x|)



o código
Manipulate[
 Plot[{
      hw[w, a],
      hw[w, 0.42]
      } // Abs // 20 Log[10, #] & // Evaluate
  , {w, -111, 111}, PlotRange -> {-120, 0}, GridLines -> Automatic, 
  PlotStyle -> {Default, Thin}, PlotPoints -> 50]
 , {{a, 0.409}, 0.3, 0.5}]








No processo de alteração do parâmetro, vamos observar algo semelhante:







Dependendo do parâmetro a, a nível do lado lobos mudanças e, a = 0,42, é mais ou menos mínima e uniforme. Com a = 0,409, podemos obter um resultado ligeiramente melhor se por “melhor” entendermos o nível mínimo possível de lobos laterais.



Nota
, , .



Desenvolvimento



Obviamente, essa microotimização não valia o esforço. É possível melhorar qualitativamente as propriedades desta janela sem alterar os dados iniciais?



Anteriormente, vimos a saída de funções de janela que somam um, o que permite restaurar o sinal original. Vamos tentar essa técnica para nossa janela e ver como ela afeta o espectro.



Vamos definir uma função auxiliar para construir janelas sobrepostas, levando em consideração os limites da mudança do argumento no intervalo de -1/2 a 1/2:







Encontre a antiderivada, desloque-a para o centro das coordenadas e dimensione-a para um:











exiba-a no gráfico:







Como você pode ver, o Wolfram fez isso sozinho aqui também, e não tivemos que especificar manualmente uma definição contínua por partes de uma antiderivada. Agora, a visualização de nossa janela depende não apenas da variável a , mas do grau de sobreposição - e à medida que aumenta, tenderá para a forma da derivada:



o código
Manipulate[
 Plot[
  OverlapShape[hsx, x, a, t]
  , {x, -1, 1}, PlotRange -> All, GridLines -> Automatic]
 , {{a, 0.404}, 0.4, 0.5}, {{t, 4}, 3, 11}]






E o toque final é encontrar uma função analítica para o espectro, a fim de determinar o valor ideal do parâmetro a .



Aqui, se tentarmos calcular a transformação diretamente, como da última vez, levaremos Wolfram a um pensamento profundo. Existem várias maneiras de acelerar este processo:



- use uma transformação discreta em vez de analítica - a mais simples. Mas um verdadeiro matemático não aplicará métodos numéricos onde uma solução puramente analítica possa ser encontrada - e não o faremos (ainda); além disso, possui limitações óbvias (em termos de resolução e limitação de espectro);



- passe valores específicos como um parâmetro de sobreposição e, olhando os resultados, tente ver o padrão e generalize-os. Uma opção bastante funcional, mas requer esforços mentais e criativos adicionais. Vamos deixar como último recurso.



- fazer cálculos diretamente no domínio da frequência. Isso é possível devido às seguintes propriedades da transformada de Fourier:



1) é linear, ou seja, a soma das imagens é igual à imagem da soma;

2) integração no domínio do tempo é equivalente a divisão por iw no domínio da frequência.

3) o deslocamento de tempo é equivalente à modulação com uma sinusóide complexa.

(mais na wikipedia ou aqui ).



Assim, tendo o espectro hw de uma função de janela arbitrária, podemos obter dele o espectro hwo de uma função somada a um com t sobreposto :











Agora podemos ver como o espectro muda na dinâmica:







Aqui, alterar o parâmetro de sobreposição já afetará o espectro da janela de uma maneira ligeiramente diferente :







Tendo brincado um pouco com os parâmetros, é fácil perceber que "mais não é melhor", e o grau ótimo de sobreposição para essa função está na região de quatro. Especificamente, para t = 4 e a= 0,404 obtemos o nível do lóbulo lateral não excedendo -80 dB. Este é um resultado muito bom - especialmente considerando que a função cosseno elevada, tradicionalmente usada para janelas somadas, fornece um nível de lóbulo de cerca de -30 dB. Bem, quando escrito explicitamente, nossa nova função de janela será semelhante a esta:



o código
OverlapShape[hsx, x, 404/1000, 4] // PiecewiseExpand // FullSimplify




{404π(12x)+36sin(13(16πx+π))+375sin(13(π8πx))606π14<x12404(2πx+π)+36sin(23π(8x+1))+375sin(13(8πx+π))606π12<x143(125cos(8πx3)+12cos(16πx3))202π+1314<x140|x|>12



Desenvolvimento adicional



O que mais você pode fazer para reduzir ainda mais o nível do lóbulo lateral? Você pode tomar os cossenos não com frequências pares, mas com frequências ímpares (aqui, para compactação, a solução do sistema de equações está embutida diretamente na definição da função):



o código
hx1 = Function[{x, a},
  a Cos[1 Pi #] + b Cos[3 Pi #] + c Cos[5 Pi #] & // #[x] /. Solve[{
          #[0] == 1,
          #[1/2] == 0,
          #'[1/2] == 0,
          #''[1/2] == 0
          }, {a, b, c}][[1]] & // # UnitBox[x] & // Evaluate]




(acos(πx)+(58a2)cos(3πx)+(38a2)cos(5πx))Π(x)



e após a integração com os parâmetros a = 0,6628 e um nível de sobreposição de 4,5 para obter uma supressão de -90 dB:







Explicitamente:

{327sin(17π(45x+2))+24855sin(17π(19x))+3670cos(114π(54x+1))+2151243024518<x1224855sin(17π(9x+1))+3670sin(37π(9x+1))+327sin(57π(9x+1))+215124302412<x5183670sin(37π(9x1))+24855sin(17π(9x+1))+327sin(57π(9x+1))43024++3670sin(37π(9x+1))+327sin(17π(45x+2))+24855sin(17π(19x))43024518<x5180|x|>12



Você pode adicionar mais um cosseno e aumentar o número de derivadas zero:



o código
hx2 = Function[{x, a},
  a Cos[1 Pi #] + b Cos[3 Pi #] + c Cos[5 Pi #] + 
       d Cos[7 Pi #] & // #[x] /. Solve[{
          #[0] == 1,
          #[1/2] == 0,
          #'[1/2] == 0,
          #''[1/2] == 0,
          #'''[1/2] == 0,
          #''''[1/2] == 0
          }, {b, c, d}][[1]] & // # UnitBox[x] & // Evaluate]




(acos(πx)+180(3516a)cos(3πx)+180(3548a)cos(5πx)+140(58a)cos(7πx))Π(x)



e após integração com os parâmetros a = 0,5862 e um nível de sobreposição de 6,4 para obter supressão de -110 dB:







Explicitamente:

{3077550sin(154π(64x5))90069sin(554π(64x5))5820sin(754π(64x5))560455cos(19π(732x))+260134452026881132<x12560455sin(118π(64x+5))+3077550sin(154π(64x+5))+90069sin(554π(64x+5))+5820sin(754π(64x+5))+2601344520268812<x11323077550sin(154π(64x5))90069sin(554π(64x5))5820sin(754π(64x5))560455cos(19π(732x))5202688++560455sin(118π(64x+5))+3077550sin(154π(64x+5))+90069sin(554π(64x+5))+5820sin(754π(64x+5))52026881132<x11320|x|>12



Uma redução ainda mais significativa no nível dos lobos laterais pode ser alcançada ao elevar a transformada de Fourier ao quadrado e, assim, reduzir o nível dos lobos laterais por um fator de 2 de uma vez. Isso permite que você se livre do aumento do número de parâmetros para sua seleção manual, mas adiciona complexidade no cálculo da convolução no domínio do tempo.



o código
hx3 = Function[{x, a}, 
  Convolve[hx1[x, a], hx1[x, a], x, y] /. y -> 2 x
 // FullSimplify // Evaluate]








e aqui já é possível atingir uma supressão acima de 160 dB:







Não daremos a fórmula explicitamente devido ao seu tamanho impressionante.



Funções de janela hipergeométrica



Para garantir o número necessário de zeros nos limites de nossa função, usamos uma busca pela solução de um sistema de equações com integração posterior. Isso não é muito conveniente - você deve alterar o número de equações todas as vezes. Talvez haja uma solução mais simples e bonita? Há sim! E a função hipergeométrica vai nos ajudar nisso.



brevemente sobre funções hipergeométricas
, . — , . — .



Nossa solução será semelhante a esta:

2zΓ(n+32)2F1(12,n2;32;z2)πΓ(n2+1)+a(z+bz3)(1z2)n/2

Onde z=sin(πx2)



Consiste na soma de duas partes: a antiderivada da função (1z2)n2normalizado em amplitude ao ponto (1,1) e uma função com parâmetros livres, multiplicados pelo mesmo (1z2)n2a fim de preservar o número necessário de derivadas zero: A







precisão do ajuste e sua influência na distribuição dos lobos laterais dependerá da função que escolhermos para corrigir a forma - há opções possíveis que requerem estudos separados. Neste caso, a funçãoa(z+bz3)oferece uma opção totalmente aceitável - e, devido a dois parâmetros, permite que você faça configurações mais precisas.



A fórmula em si é interessante (e conveniente) porque para n par é simplificada para um polinômio e para n ímpar é simplificada para a soma de um polinômio com arco-seno e raiz quadrada - o que torna a fórmula final mais compacta e fácil de calcular.



Já é mais fácil (e rápido) selecionar os parâmetros aqui por meio da transformada discreta de Fourier. Também precisamos de uma definição adicional da função de sobreposição para trabalhar com três parâmetros:



















Por exemplo, após substituir os parâmetros da imagem, nossa função será simplificada para



o código


288z11+2315z97380z7+12330z511940z3+8163z3200





Janela Kaiser Inversa



Se a tarefa de somar as janelas em uma não valer a pena, então a janela ideal é a janela Kaiser, que permite ajustar suavemente a largura do lóbulo principal. No entanto, tem uma desvantagem - uma vez que é expresso em termos da função de Bessel , é um tanto difícil lê-lo fora dos pacotes matemáticos. Você pode, é claro, implementar separadamente a função de Bessel - ou pode procurar uma aproximação por meio de funções elementares. E de repente descobriu-se que usando a função do espectro da janela Kaiser (limitada no tempo) para isso, você pode obter um resultado ainda um pouco melhor -



sinh(k14x2)sinh(k)14x2Π(x)

Nota
±12 , ksinh(k).


O pagamento por uma decisão tão astuta foi a impossibilidade de expressar seu espectro de forma puramente analítica - mas isso não é tão importante, já que na prática, em qualquer caso, operamos e analisamos dados de forma discreta usando a transformada discreta de Fourier. Abaixo no gráfico você pode comparar seus espectros, onde vermelho é a janela original do Kaiser, verde é sua aproximação:



o código








Com a mesma largura da pétala principal, o nível dos lobos laterais é ligeiramente inferior. E para facilidade de uso, você pode elaborar uma tabela com valores aproximados do parâmetro k para obter o nível necessário de supressão do lóbulo lateral:

-60 k=8.8-90 k=11.36-120 k=15.18-150 k=18.88



Encontrar uma aproximação da antiderivada da janela Kaiser inversa permanece um problema interessante separado. Até agora, só foi possível expressá-lo por meio de uma série de potências, que convergem um tanto lentamente:



sinh(k14x2)sinh(k)14x2dx=n=1212nkn1x2n1(k(1)nKn12(k)+kKn12(k))π(2n1)sinh(k)Γ(n)



Talvez os especialistas em matemática possam sugerir uma solução melhor nos comentários.



Conclusão



Apesar do processamento digital de sinais ser uma disciplina completamente formada, ainda há espaço para desenvolvimento e algo novo nele - pelo menos na literatura científica e especial existente, as questões de construção de funções de janela resumidas em uma unidade não são consideradas; e as capacidades dos modernos sistemas de álgebra computacional permitem repensar o conhecimento acumulado em um nível qualitativamente novo.



O documento original do Wolfram Mathematica com todos os cálculos está disponível no GitHub .



All Articles