ANÁLISE ESPECTRAL - CONTINUACÃO

ROTEIRO PARA ANÁLISE ESPECTRAL

A discussão a seguir exemplifica um possível roteiro para o cálculo da análise espectral. Este roteiro pode ser acompanhado juntamente com o programa spectrum.good.f fornecido para este propósito. Notem que o programa esta preparado para fazer analise espectral em um período de 120 dias por 21 anos. Você pode modificar o programa de acordo com as suas necessidades. O programa foi desenvolvido por Charles Jones (cjones@icess.ucsb.edu).

 Este método utiliza Fast Fourier Transform (FFT) para o seu cálculo e pode ser compilado e executado, desde que o computador tenha o pacote starpac que esta disponível no NCAR em http://www.scd.ucar.edu/softlib/mathlib.html

O nome da rotina do starpac utilizada no programa e' PGMS

 Pode-se obter o espectro de outras formas, como por exemplo, utilizando-se a função de auto-covariância (ex. Mitchell, 1966).  Na presente discussão, iremos nos concentrar no roteiro que utiliza FFT em seu procedimento.

Explicações sobre as variáveis utilizadas estão comentadas dentro do próprio programa, na maioria das vezes. Para mais detalhes sobre as variáveis , consultem:

Starpac.guide.doc

Usem <editar> <procurar> ou <edit> <search> para achar o significado de variáveis de entrada não comentadas no programa (caso julguem necessário)

********************************************************************************************************************************

OBS: Notem no programa que uma das variáveis de entrada e' o NFFT. Abaixo vai uma explicação sobre a mesma extraída do starpac_guide.doc

NFFT --- The minimum extended series length that meets the requirements of
the Singleton FFT code. NFFT is returned by FFTLEN; it must be
input to PGMS and FFTR; it is computed internally by all of the
other subroutines using the FFT.

The extended length of the series is

NFFT = N + LAGMAX + K

where N is the number of observations in the series. K >= 2 is
chosen so that NFFT is as small as possible, NFFT-2 is divisible
by 4 and has no prime factors greater than 23, and the product of
the square-free prime factors of NFFT-2 does not exceed 209. In
general, NFFT will be less than N+LAGMAX+100 for the correlation,
Fourier spectrum and autoregressive spectrum subroutines. STARPAC
subroutine FFTLEN, when called using the sequence
 

LAGMAX  --> The maximum lag  value for which the correlation  coefficients are
             computed.   The default value of LAGMAX is selected  by STARPAC as
             follows.

               LAGMAX =  100      if   301 <= N.

               LAGMAX =  N/3      if    96 <= N <= 300.

               LAGMAX =   32      if    33 <= N <=  95.

               LAGMAX =  N-1      if          N <=  32.

             If LAGMAX is less than or equal  to zero or if neither  LAGMAX nor
             LAGS is an argument in  the CALL  statement the  default  value is
             used.   When LAGS is an argument in the  CALL statement, LAGMAX is
             the largest value in LAGS.
 

******************************************************************************************************************************************

EXEMPLO DE UTILIZAÇÃO DO ROTEIRO PARA A OBTENÇÃO DE ESPECTRO

Considere uma série temporal y = y(t), t = 1,N. A série temporal pode ser, por exemplo, medias diárias da temperatura do ar próximo à superfície na cidade de São Paulo. Os valores diários vão de 1 de janeiro de 1979 a 31 de dezembro de 1979 (portanto N=365).

1)     O primeiro passo consiste em remover possíveis tendências lineares na serie temporal y(t). Portanto ajusta-se uma reta z = a + bt a série temporal y(t). Em seguida efetua-se a operação:

                                          (1)      

2)     A nova série y’(t) pode conter flutuações com períodos já conhecidos. No exemplo em questão, a série original y(t) pode exibir um ciclo anual, que também estará presente na nova serie y’(t). Se o objetivo em questão e’ estudar o espectro em freqüências próximas ao ciclo anual, recomenda-se remover o ciclo anual da série y’(t). Este passo consiste então realizar uma analise harmônica na serie y’(t) e reter o primeiro e segundo harmônicos:

                                           (2)

A função f(t) é subtraída da série temporal y’(t): y”(t) = y’(t) – f(t).

Se  e’ a média de y”(t), em seguida remove-se  de y”(t): y’’’(t) = y”(t) -  .

• OBS: O primeiro e o segundo Harmônicos só representam o ciclo anual e o semi-anual somente porque temos uma série temporal de 365 dias. Lembrem-se que quando k=1 (primeiro harmônico) estamos assumindo que há um ciclo completo no intervalo de tempo total da série (no caso N=365 dias). Se a série for maior do que o comprimento de um ano, não será mais o primeiro harmônico o que representa o ciclo anual.

3)     O passo seguinte consiste em aplicar um “taper” nas extremidades da série temporal y’’’(t). O procedimento de “tapering” é recomendado para evitar uma descontinuidade no final da série de dados. Esta descontinuidade poderia representar a adição de um cosseno que na realidade não existe. Em geral, utiliza-se o taper em 5% (ou 10%) em cada extremidade. Varias tipos de funções taper são disponíveis. Uma função taper bastante utilizada é chamada de “split-cosine-bell”. A operação consiste em:

y*(t) = w(t) * y’’’(t)                                          (3)

onde a função w(t), split-cosine-bell, e’ dada por:

w(t) = 0.5*(1-cos[pi*(t-0.5)/m])                                1 <= t <= m

w(t) =   1                                                             m + 1 <= t <= N-m (4)

w(t) = 0.5*(1-cos[pi*(N-t+0.5)/m])                N-m + 1 <= t <= N

e m é calculado tal que 2m/N = TAPERP é a proporção da série a ser aplicado o taper (e.g. 5% ou 10%). Ver discussão em Bloomfield (1976) para maiores detalhes da importância de se aplicar o taper antes de calcular a analise espectral.

4)     A série y*(t) resultante do passo anterior e’ utilizada para se calcular o espectro. Aplica-se a transformada de Fourier rápida (FFT) para o cálculo das estimativas “brutas” das amplitudes espectrais (“raw spectral estimates”). Os resultados deste passo são as amplitudes espectrais C²_k, k = 1,N/2.

5)     As estimativas “brutas” das amplitudes espectrais são estatísticas enviesadas (com bias) do espectro. Isto e, se mais pontos forem adicionados ou removidos da serie, o espectro pode ser alterado. Dessa forma, o passo seguinte consiste em alisar as estimativas espectrais. Varias janelas espectrais são disponíveis (Ver discussão em Bloomfield 1976 e Chatfield 1996). Alguns exemplos de janelas são: Turkey, Parzen, Hamming, Daniell (veja discussao no livro do Chatfield, pag 116, ed 1996). Um método bastante comum consiste em alisar o espectro com uma media móvel de comprimento L com igual peso para cada ponto (método Daniell). Em seguida plota-se o espectro de freqüência.

• OBS: A escolha de L depende das características do espectro as quais dependem do tanto de freqüências que se esta resolvendo. Por exemplo, quando temos uma serie longa de, digamos, 20 anos de dados diários, temos muitas freqüências a se resolver. Neste caso, quando calcularmos o espectro, veremos que o mesmo tem um aspecto muito 'ruidoso'. Para eliminar estes ruídos temos que utilizar um L maior (de preferência maior do que 3). Entretanto, podemos estar interessados em obter o espectro de uma serie temporal de dados diários restritos apenas aos meses de verão. Uma forma de se obter este espectro e' obter o espectro a cada ano e então obter a media do mesmo. Neste caso, ao calcularmos a media já estamos 'alisando' o espectro e L pode ser menor (por exemplo, igual a 3). De qualquer maneira, o bom-senso deve sempre imperar e um bom julgamento de quem esta efetuando a analise e' recomendado.

6)     Estimar o número de graus de liberdade. Este cálculo necessita levar em consideração a aplicação do taper antes da análise espectral. Para uma discussão mais detalhada ver Madden and Julian (1971) ou Madden (1977). De acordo com os autores, uma melhor estimativa dos graus de liberdade quando se procede ao taper, e' obtida definindo-se um comprimento equivalente da serie Neff dado por:

     Onde N e' o comprimento da serie (estamos considerando aqui o taper em 10% no inicio e no final da serie)

     O numero de graus de liberdade e’ dado aproximadamente por:

Dof @ 2* L * Neff/N = 2 *L * 0.873                                        (5)

Notem, entretanto, que se o espectro for calculado como uma media (ensemble) de NS espectros, então Dof deve ser aumentando, ficando:

Dof @ 2* L * (Neff/N)*NS = 2 *L * 0.873 * NS                                        (6)

7)     O passo final consiste em avaliar o espectro de fundo e determinar se o espectro calculado difere do espectro de fundo por uma quantia estatisticamente significante. Começa-se primeiramente em ajustar uma “hipótese nula” para o espectro. Em geral, se auto-correlação de lag 1 (R₁) não difere de zero, a série temporal y*(t) é considerada como livre de persistência. Neste caso, uma estimativa apropriada para a hipótese nula é que o espectro de fundo e’ dado por um “ruído branco”. Em outras palavras, o espectro de fundo consiste em uma linha paralela ao eixo de freqüências e interceptando o eixo das abscissas em , onde  é a média das amplitudes espectrais.

Se por outro lado, R₁ difere de zero e as auto-correlações seguintes seguem uma aproximação exponencial: R₂ @ R₁², R₃ @ R₁³ etc então o espectro de fundo é denominado de “ruído vermelho” e dado por:

   (6)

O espectro de fundo (hipótese nula) e’ calculado e superposto no gráfico do espectro calculado no passo (5) acima.  O nome 'ruído vermelho' ou 'ruído branco' vem da analogia com o espectro da radiação. Quando todos os comprimentos de onda tem a mesma importância, temos a luz branca como resultado. Quando os comprimentos de onda mais longos (menores freqüências) possuem maior amplitude, o espectro move-se para o vermelho.  Na atmosfera, e' muito comum a ocorrência do espectro de fundo ser do tipo 'ruído vermelho' devido as características de 'curta memória' dos processos físicos.

8) Em seguida, deseja-se determinar se as amplitudes espectrais diferem das amplitudes do espectro de fundo. Tukey (1950) determinou que a razão entre a amplitude espectral do espectro calculado e amplitude do espectro de fundo segue uma distribuição estatística c² dividida pelo número de graus de liberdade (dof):

                                     (7)

Em realidade, a Eq. (7) indica uma igualdade de distribuições. Na prática, dado o número de graus de liberdade e o nível de significância desejado (por exemplo, 5%), usa-se a tabela de distribuição de freqüência c² para se determinar o valor de corte. Este valor de corte deve ser dividido pelo numero de graus de liberdade. Para obtermos a curva do 95% de confiança, basta multiplicarmos as ordenadas do espectro de ruído vermelho pela razão entre o valor de corte e o DOF. Por exemplo, se DOF=30 para um intervalo de confiança de 95% (equivalente a nível de significância de 5%) o valor de corte e' 43.773/30 = 1.459 (43.773 e' extraído da tabela c²). Clique no link abaixo para poder observar a tabela da probabilidade c²/ν  (ν=Dof). Para utilizá-la basta encontrar Dof que ja se obtém o lado esquerdo da Eq. 7

TABELA c²/ν

  Um exemplo de aplicação esta sendo mostrado na figura abaixo abaixo. No painel da esquerda mostra-se o vento médio e a vorticidade relativa em 700hPa obtidos de Junho a Setembro (1979-2000) com dados diários. Na figura do vento médio estão destacadas 3 caixinhas A, B, C. Em sentido horário na Fig. 1 podemos observar os espectros calculados para cada uma das caixas. As curvas do ruído vermelho (linha cheia) e de 95% de significância (linha tracejada) estão também mostradas. Conclui-se destas figuras que apenas as amplitudes espectrais que superam a linha tracejada são significativas ao nível de significância de 95%. Os demais picos são resultado do fato da atmosfera apresentar 'memória' e são considerados parte do ruído-vermelho que decorre deste fato.  Para maiores detalhes, ver discussão em Mitchell (1966).

Figure 1. Top left : seasonal mean wind field at 700 hPa. Wind speed is indicated with shading (m s^-1). Seasonal mean winds at 700 hPa (vectors) and relative vorticity at 700 hPa (contours). Solid (dashed) contours denote positive (negative) relative vorticity with intervals of 0.5 x 10^-5 s^-1. Summer season is defined as 1 June–30 September 1979-2000. Boxes indicated with thin solid lines (A, B, C) are used to compute spectral variance. Other panels indicate spectral variance of relative vorticity at 700 hPa at regions A (top), B (middle) and C (bottom) indicated by thin solid lines on top left (From Jones et al. 2003).

• LARGURA DE BANDA (bandwidth):

    O espectro final obtido foi alisado com um dos procedimentos descritos acima (por exemplo, Daniell - media movel com peso igual para todos o pontos). Quando o periodograma e alisado por uma media movel de tamanho L, esta media produz uma estimativa de um espectro continuo visto atraves de uma JANELA RETANGULAR de largura de banda igual a (L/(N/2))*fn onde L e' o comprimento da media movel, N o numero de pontos (ou N/2 o numero de frequencias resolviveis possivel) e fn e' a frequencia de Nyquist. Assim, o calculo da largura da banda (bandwidth) e' importante para que possamos identificar quao realista e' a separacao de 2 picos significativos em um espectro. A largura de banda pode ser dada tanto em frequencia quanto em periodo. Assim, dada a bandwidth, sabemos que 2 picos que se encontrem em distancia inferior a essa largura nao sao, de fato, separaveis, dado o processo de de alisamento do espectro (veja Madden 1977 para detalhes).

******************************************************************************************************************

ROTEIRO PARA CÁLCULO DE ESPECTRO UTILIZANDO O STATISTICA 6.0 (EXERCÍCIO SOBRE ESPECTRO)

      1)      Salve os dados que deseja utilizar numa planilha excel. Lembre-se de deixar apenas a primeira linha com o nome das variáveis (no caso dos dados de TSM, as coordenadas)

2)       Abra o Statistica 6.0, clique em

• File

• Open (escolha opção xls no caso de dados excel. Abra, por exemplo, os dados de TSM disponíveis no site)

• Import select sheet  (confira o número de linhas e colunas)

• Get variable names from first row (clicar nessa opção para pegar o nome das variáveis na primeira linha)

 3)       Com a planilha aberta selecione a coluna que deseja calcular o espectro clicando no topo da mesma. Note que a coluna inteira ficará marcada

4)      Vá no menu do topo da página e clique em:

·        Statistics

·        Advanced Linear/Nonlinear modes

·        Time series/Forecasting

    5)      Aparecerá, então uma janela com os dados que você escolheu mais diversas opções de análise. Clique em:

·        Spectral (Fourier) analysis

6)       Ao clicar na opção acima você entrará em uma janela que contém diversas opções. LEMBRE-SE QUE SEMPRE QUE VOCÊ TIVER DÚVIDAS SOBRE O QUE ESTÁ NA JANELA, CLIQUE EM ? NO TOPO DA JANELA QUE UM HELP SERÁ ABERTO EXPLICANDO AS OPÇÕES (OU PELO MENOS A MAIORIA DELAS).

Principais opções a serem consideradas (veja também explicação da matéria em aula):

·        Taper (clique na opção e escolha 10% - e.g. Madden (1977))

·        Subtract mean (ative)

·        Detrend (ative)

·        Pad length to power of 2 (e.g. JW, CH) (importante para o procedimento de FFT)

·        Single series Fourier analysis ( finalmente, clique nessa opção)

7)       Uma nova janela deve aparecer com algumas opções. O principal agora é decidir a janela para a estimativa da densidade espectral. Você pode experimentar diferentes janelas com diferentes métodos e checar seus resultados. As principais teorias sobre as janelas encontram-se, por exemplo, em CH, cap. 7. No presente caso, vamos adotar a média móvel simples, que considera pesos iguais para todos os pontos da janela. Um exemplo dessa utilização está em Madden (1977) e Madden (1971). A opção de pesos iguais para todos os pontos é a opção Daniell

·        Width of data window: (essa decisão depende um pouco do caso, mas vamos começar escolhendo 5 pontos – número de pontos da média-móvel.- e.g. Madden 1977)

·        N largest values (escolha arbitrária. Na realidade, a decisão será feita numa análise posterior de verificação de significância)

·        Density estimates (clique nessa opção). Demais opções, ficam a critério de cada um

·        Summary (clique nessa opção para ver um sumário dos resultados)

Para retornar para as opções iniciais, é só clicar novamente na barrinha “Single Series Fourier” que se encontra na barra inferior da janela maior do statistica.

8)       Plote os resultados e reflita sobre os mesmos. Tente as diferentes opções que aparecem (como plotar com o período, com a freqüência, plotar a densidade espectral ou o periodograma).

9) O periodograma e' definido como a soma do quadrado dos coeficientes do seno e do co-seno multiplicado por n/2. A densidade espectral e' calculada apos aplicar a janela espectral (no caso de Daniel, sera uma media movel de comprimento K)

10)     Salve os resultados e faça a análise de significância conforme visto no roteiro da sala de aula, agora usando o excel e uma tabela do χ2. Lembrem-se ainda de calcular os graus de liberdade, conforme visto no roteiro. Voce pode utilizar o periodograma ou a densidade espectral (ou ambos) para comparar seus resultados. Para salvar, simplesmente faça o seguinte: Com a janela ativa com a tabela dos resultados, clique em

• Edit

• Copy with headers

Abra uma planilha excel nova e faça o “paste” com <coltrol><c>

10)  Determine também a Bandwidth (no caso da média móvel com iguais pesos – ou o procedimento “Daniell”, é dada por :

            Onde L é o comprimento da média móvel (5 no caso), N é o comprimento da série (ou N/2 o número de freqüências possíveis de se resolver) e F_Nq é a freqüência de Nyquist

Nota: Para determinar o espectro de fundo, você precisará calcular a auto-correlação para diferentes lags. Você pode fazer isso fazendo os seguintes passos:

i)                    Usando o Statistica, Retorne ao nível anterior ao cálculo do espectro.

ii)                  Escolha a opção  ARIMA & AUTOCORRELATION FUNCTIONS

iii)                Clique na opção AUTOCORRELATION (com WHITE NOISE acionado)

iv)                Veja se o correlograma difere do white noise

v)                  Salve os dados e decida se a aproximação de ruído-vermelho para o espectro de fundo é válida

FIM DOS PROCEDIMENTOS

Referências

Bloomfield, P.: 1976: Fourier analyses of time series: an introduction. New York.

Chatfield C., 1996: The Analysis of Time Series: An introduction. Chapman & Hall, fifth edition, NY. 283 pp.

Jones, C. N. Mahowald, and C. Luo, 2003: The role of easterly waves on African desert dust transport. J. Climate (in press).

Madden, R. A., and P. R. Julian, 1971: Detection of a 40-50 day oscillation in the zonal wind in the tropical Pacific. J. Atmos. Sci., 28, 702-708.

Mitchell, J. M. jr., 1966: Climate Change. World Meteorological Organization, Geneva, Tech. Note no. 79, 79 pp.

Roteiro analise espectral .doc           Roteiro para usar o estatistica .doc

[ Home ] [ Up ]

Send mail to leila@model.iag.usp.br with questions or comments about this web site.
Last modified: 11/08/05