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:
- 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 ...
- 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 ...
- 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:
- sem lacunas nas bordas - valores iguais a zero. As pausas são de Dirichlet, Hamming, Blackman-Harris. Não tem - Triangular, Hann, Nuttal
- ausência de descontinuidades do 1º e outros derivados
Mais duas funções de janela se destacam separadamente:
- Kaiser, que permite definir a largura da pétala principal;
- 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
↓
NuttallWindow[x] // FunctionExpand
↓
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}]
↓

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:
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
↓
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]
↓
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:
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]
↓
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:
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:
Onde
Consiste na soma de duas partes: a antiderivada da função normalizado em amplitude ao ponto (1,1) e uma função com parâmetros livres, multiplicados pelo mesmo a 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çãooferece 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
↓

↓
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 -
Nota
, .
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:
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:
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 .