KR100817692B1 - 이산 푸리에 변환에 의한 시계열 데이터 위상 추정 방법 - Google Patents

Abstract

본 발명은 푸리에 변환에 의한 시계열 데이터의 위상 추정 방법에 관한 것으로, 푸리에 스펙트럼이나 전력스펙트럼에서 첨두치와 첨두 주파수 및 모드의 위상 사이에 일정한 수학적 관계에서 신호의 위상을 추정할 수 있다. 시계열 데이터에 이산 푸리에 변환을 적용한 결과 얻어진 푸리에 스펙트럼에서 첨두치를 중심으로 좌우 스펙트럼을 복소지수감쇠함수로 적합함으로써 모드를 추정할 수 있고, 각 모드에 대응하는 위상을 추정할 수 있다. 본 발명은 단순 산술 계산 과정과 첨두치 탐색 과정으로 신호 모드에 대응하는 위상을 추정할 수 있으므로, 계산시간이 매우 빠르고 실시간 대용량 신호처리를 가능하게 한다.

Description

이산 푸리에 변환에 의한 시계열 데이터 위상 추정 방법 {A method for estimating phase angle of time series data by discrete Fourier transform}

도 1은 본 발명을 수행하기 위한 하드웨어의 개략적인 블록 구성도,

도 2는 본 발명의 실시예의 동작 과정을 간략히 도시한 플로우챠트,

도 3은 비주기 신호의 푸리에 스펙트럼 크기,

도 4는 이산신호에 대한 첨두주파수와 푸리에 스펙트럼 크기,

도 5는 이산신호에 대한 첨두주파수와 푸리에 스펙트럼 위상,

도 6은 실시예에 대한 이산 시계열 신호,

도 7은 실시예에 대한 푸리에 스펙트럼 크기,

도 8은 실시예에 대한 푸리에 스펙트럼 위상,

도 9는 이산 푸리에 변환 결과 스펙트럼이 비대칭인 경우,

도 10은 이산 푸리에 변환 결과 스펙트럼이 대칭인 경우.

<도면의 주요 부분에 대한 부호의 설명>

10: 입력 처리부 20: 연산부

30: 출력부 40: 이산시간 발생부

본 발명은 푸리에 변환에 의한 시계열 데이터의 위상 추정 방법에 관한 것으로, 이산 시계열 데이터에 포함된 모드나 파라미터를 추정하는 시스템 동정(system identification) 이나 통계적 신호처리(statistical signal processing)와 관련되어 있다.

지금까지 시계열 데이터에서 파라미터를 추정하는 방법은 크게 파라메트릭 (parametric) 방법과 비파라메트릭(non-parametric) 방법으로 구분할 수 있다. 파라메트릭 방법은 신호를 자기회귀이동평균(autoregressive moving average, ARMA) 모델로 변환하여 얻은 선형예측행렬(linear prediction matrix)로 부터 복소지수함수를 계산하여 이산 데이터에 포함된 파라미터(모드)를 추정하는 방법이다.

대표적인 파라메트릭 방법으로 프로니 법(Prony method)이 있는데, 이 알고리즘은 이산신호를 복소지수함수로 적합할 때, 선형예측방정식을 풀어야 하고, 고차방정식의 해를 계산해야하므로 계산 시간이 많이 소요되며 대용량 기억장소를 가진 컴퓨터를 필요로 한다. 또한 이 방법은 잡음에 대해서 오차를 수반하며 신호데이터의 구간간격과 데이터의 크기에 따라서 파라미터가 민감하게 변화하는 단점이 있다.

이산 푸리에 변환은 다양한 산업분야에 널리 적용되고 있고 주로 신호에 포함된 주파수를 계산하는데 이용하고 있다. 이와 같은 이산 푸리에 변환을 이용하는 비파라메트릭 방법은 푸리에 변환을 반복적으로 적용하여 이산 데이터에 포함된 파라미터를 추정하는 방법이다. 비파라메트릭 방법은 반복적으로 이산 푸리에 변환을 적용하여 각 주파수에 대응하는 진동모드를 추정하는 방법으로 각 첨두치에 대응하는 스펙트럼의 크기를 연속적으로 계산한 후 곡선적합을 적용하여 첨두주파수에 대응하는 모드를 추정한다. 그러나 몇 번의 푸리에 변환을 수행해야 하고, 스펙트럼 누설이 발생하면 오차가 크게 발생하는 단점이 있으며, 특히 본 발명에서 구현하고자하는 각 모드에 대응하는 위상을 추정하는 방법은 지금까지 개발된 바 없다.

상기한 종래의 방법은 선형방정식, 최소자승법과 같은 복잡한 수치적인 알고리즘을 포함하고 있어 계산시간이 많이 소요될 뿐만 아니라 큰 컴퓨터 기억장소를 필요로 한다.

이에 본 발명은 종래기술의 문제점을 해소하기 위해 이산 푸리에 스펙트럼으로부터 직접 모드의 위상을 추정할 수 있는 방법을 제공하는 데 있어, 이산 푸리에 변환 알고리즘을 이용하므로 신뢰성이 있고, 계산 과정이 매우 단순하므로 효율적이며, 계산 시간을 크게 단축할 수 있을 뿐만 아니라 소용량의 기억장소를 가진 컴퓨터로 시계열 데이터에 포함된 각 모드에 대한 위상을 추정할 수 있는 방법을 제공하는 데 그 목적이 있다.

[본 발명을 수행하기 위한 하드웨어]

도 1은 본 발명을 수행하기 위한 하드웨어의 개략적인 블록 구성도로서, 동도면을 참조하면 알 수 있듯이, 본 발명을 수행하기 위한 하드웨어는, 입력 처리부(10)와 연산부(20)와 출력부(30)와 이산시간 발생부(40)를 포함하여 구성될 수 있다.

상기 입력 처리부(10)는 입력되는 아날로그형태의 시계열 데이터(time series data)를 아날로그 투 디지털 변환(analog to digital convert)하여 연산부(20)로 입력한다.

상기 연산부(20)는 입력 처리부(10)로부터 입력되는 시계열 데이터에 대해 이산 푸리에 변환(discrete Fourier transform)을 수행하여 모드(mode)를 추정하고 이로부터 각 모드에 대응하는 시계열 데이터의 위상을 추정한다.

상기 출력부(30)는 상기 연산부(20)에 의해 추정된 위상의 추정 결과를 출력한다.

상기 이산시간 발생부(40)는 상기 연산부의 푸리에 변환을 위해 필요한 시간 정보인 이산시간 데이터를 발생하여 이를 연산부(20)로 입력한다.

이하에서는, 상기한 바와 같이 구성된 하드웨어를 통해 수행되는 본 발명의 실시예에 대해 설명하기로 한다.

본 발명은, 시계열 데이터의 이산 푸리에 변환(discrete Fourier transform)에서 시스템의 상태를 나타내는 각 모드에 대한 위상을 직접 구할 수 있는 방법에 관한 것으로, 푸리에 스펙트럼으로부터 위상을 추정하는 실시예를 포함한다.

우선, 본 발명의 실시예에 적용되는 원리에 대해 설명하기로 한다. 연속함수에 대한 푸리에 변환 식은 수학식 1과 같다.

실제 시계열 데이터로 주어진 신호는 주기적인 신호보다 비주기 신호인 경우가 대부분이다. 구간 0≤n≤N-1에 정의된 신호 x[n]의 이산 푸리에 변환은 수학식 2와 같다.

수학식 2에서 W_N =exp(-j2π/N)이며, 스펙트럼 계수 X[k]는 구간 0≤n≤N-1에서 정의된다. 한편 신호가 비주기적이라면 이산 푸리에 변환의 스펙트럼은 연속적으로 표현되는데, 대부분 자연계의 신호는 비주기적이므로 스펙트럼은 연속적으로 표현된다. 따라서 이산 푸리에 변환 결과 얻어진 푸리에 스펙트럼에서 첨두 스펙트럼(peak spectrum)을 중심으로 한 좌우 스펙트럼은 포락선 형태가 된다[도 3].

따라서 본 발명에서는 푸리에 스펙트럼에서 첨두 스펙트럼을 포함하여 좌우스펙트럼을 복소지수감쇠함수(수학식 3 또는 수학식 4)로 가정하였다. 그러므로 복소지수감쇠함수의 푸리에 스펙트럼으로부터 신호에 포함된 주파수와 제동계수를 추정할 수 있고 각 모드에 대응하는 위상을 추정할 수 있다[도 2].

연속함수 x(t)가 주파수 ω₁에서 제동계수 α와 크기 A, 위상 Φ₁를 가진 지수감쇠코사인함수라 하면 수학식 3과 같이 나타낼 수 있다. 그리고 위상 Φ₂를 가 진 지수감쇠사인함수라 하면, 수학식 4와 같이 나타낼 수 있다.

먼저 지수감쇠코사인함수에 대해서 고려하기 위해서 수학식 3을 푸리에 변환하면 수학식 5와 같다.

주파수 ω₁에 대응하는 위상을 Φ₁이라 하고, 제동계수가 작아서 ω₁≫α이 성립할 때, 수학식 6과 같다.

이 식을 정리하면, 수학식 7과 같이 간단하게 표현된다.

그러므로 스펙트럼의 크기 X_ω1과 각도∠X(ω₁)을 분리하여 표현하면, 수학식 8, 수학식 9와 같다.

따라서 주파수 ω₁에 대응하는 위상 φ₁을 가진 지수감쇠코사인함수에서 위상은 푸리에 스펙트럼의 위상과 같다.

한편 수학식 4에 나타난 지수감쇠사인함수를 코사인함수로 표현하면 다음과 같다.

그러므로 수학식 9로 표현된 지수감쇠코사인함수의 위상과 스펙트럼 위상 관계에서 지수감쇠사인함수의 위상은 수학식 11과 같이 표현할 수 있다.

따라서 지수감쇠사인함수의 위상과 푸리에 스펙트럼 위상 사이에는 수학식 12와 같은 관계가 성립한다.

수학식 12에서 지수감쇠사인함수의 위상은 푸리에 스펙트럼의 위상보다 위상이 90°지연됨을 알 수 있다.

그러므로 이산 시계열 데이터에 포함된 파라미터에서 위상은 주어진 시계열 데이터에 대한 이산 푸리에 변환으로부터 직접 추정할 수 있다. 임의의 이산 시계열 데이터를 이산 푸리에 변환하면, 도 4와 같이 푸리에 스펙트럼의 크기와 도 5와 같이 푸리에 스펙트럼 위상으로 나타낼 수 있다. 그리고 첨두 주파수 ω₁이 중요 모드에 대응하는 주파수라면, 위상은 φ₁이 되어 매우 쉽게 위상을 추정할 수 있다.

다음은 실시예를 통하여 본 발명의 정확성과 효율성을 검증한다. 도 6은 시간함수 y(t)=1.0e^{-0.1 t}cos(4.2t+30°)에 대해서, 샘플링 구간을 1/60초로 설정하고, 30초 동안 샘플링 한 데이터(N=1796)를 나타낸 것이다. 그림에서 지수적으로 감쇠하고 있는 정현파임을 알 수 있다. 그리고 도 7은 도 6에 나타난 신호에 이산 푸리에 변환을 적용한 결과 얻은 전력스펙트럼을 도식한 그림으로 푸리에 스펙트럼의 크기와 유사한 형태를 가지고 있다. 이로부터 첨두주파수가 0.67Hz(4.2[rad/sec]) 근처에 존재함을 알 수 있다. 즉, 전력 스펙트럼의 첨두치가 0.67Hz 근처에 형성되어 있어 중요 모드의 주파수가 0.67Hz임을 알 수 있다.

또한 도 8에는 푸리에 스펙트럼의 위상을 도식한 것이다. 첨두 스펙트럼에 대응하는 주파수에 대한 위상은 약 40°아래에 형성되어 있고, 수치적 계산에 의하면 정확하게 32.2°를 추정하고 있어 초기 주어진 위상 30°에 근접한 결과를 산출하였음을 알 수 있다. 즉, 중요 모드에 대응하는 위상은 스펙트럼의 위상과 같으므로 도 8에 나타낸 것과 같이 단순한 복소수 계산에서 중요 모드의 위상을 추정할 수 있음을 알 수 있고 본 발명의 효율성을 확인할 수 있다.

한편, 이산 푸리에 변환은 주어진 데이터의 수가 적을 때, 스펙트럼의 누설이 발생한다. 이산 푸리에 변환 결과, 첨두치 주변에 영향을 크게 끼치는 모드가 없다고 가정할 때, 스펙트럼 누설이 크게 발생하면, 첨두치 좌우 스펙트럼은 첨두값을 중심으로 비대칭이 크게 되고[도 9], 누설이 작으면, 첨두치를 중심으로 좌우 스펙트럼은 대칭이 된다[도 10]. 즉, 도 9와 같이 첨두스펙트럼의 좌우 스펙트럼이 비대칭일 때, 첨두주파수가 정확한 값이 아니므로 첨두주파수에 대응하는 위상도 정확한 값이 아니다. 그러나 도 10과 같이 첨두스펙트럼의 좌우스펙트럼이 대칭에 가까우면, 첨두주파수는 정확한 값이므로 이를 이용해서 계산한 위상은 정확한 값이 된다. 따라서 푸리에 스펙트럼 위상으로부터 추정한 모드의 위상의 정확성 판단은 스펙트럼의 대칭성으로부터 판단할 수 있다. 도 9에서 첨두치 좌측 스펙트럼(S210)의 크기 X₁과 우측 스펙트럼(S220)의 크기 X₂를 이용하여 대칭 정도를 정량화하면 수학식 13과 같다. 스펙트럼의 크기가 X₂ ＞ X₁ 일 때,

만일 수학식 14와 같이 새로운 변수를 정의하면,

수학식 13은 수학식 15와 같이 스펙트럼 편차의 함수로 정리된다.

수학식 15에서 ΔX_idx가 0%에 가까울수록 스펙트럼으로부터 추정한 위상은 정확한 값이라 판정할 수 있고, ΔX_idx가 100％에 가까울수록 추정한 위상은 부정확한 값이라 판정할 수 있다.

상술한 바와 같이 본 발명은, 시계열 데이터의 푸리에 변환에서 직접 동적시스템의 특성을 나타내는 위상을 추정할 수 있는 것으로, 기존의 방법들과는 달리 반복적인 시뮬레이션에 의존하지 않기 때문에, 푸리에 변환을 반복 계산하지 않고, 1회의 이산 푸리에 변환에서 신호에 포함된 각 모드에 대응하는 위상을 추정할 수 있다.

즉, 본 발명은 위상 추정 시간을 크게 단축하므로, 실시간(real time)으로 빠른 파라미터 추정을 요구하는 시스템에서 성능 향상 등을 이룰 수 있다. 특히, 본 발명은 수학식에 근거를 둔 단순하면서 명확한 알고리즘이므로 이산 푸리에 변환을 적용하는 모든 시스템에 용이하게 적용할 수 있는 효과가 있다.

Claims (4)

1. 이산 푸리에 변환에 의한 시계열 데이터의 위상 추정방법에 있어서,

   데이터취득 장치로부터 이산치 시계열 데이터를 취득하는 스텝과,

   상기 취득된 이산치 시계열 데이터에 대해서 이산푸리에변환을 수행하여 이산치 푸리에 스펙트럼을 계산하는 스텝과,

   상기 계산된 이산치 푸리에 스펙트럼으로부터 첨두치를 선택하고, 첨두치에 대응하는 푸리에 스펙트럼 위상을 계산하는 스텝과,

   상기 계산된 푸리에스펙트럼 위상을 시계열 데이터에 포함된 첨두주파수에 대한 위상으로 치환하는 스텝과,

   이 치환 결과를 출력하는 스텝을 포함하는 것을 특징으로 하는 이산푸리에 변환에 의한 시계열 데이터의 위상 추정방법.
2. 청구항 1에 있어서,

   상기 이산치 푸리에 스펙트럼의 위상을 시계열 데이터에 포함된 첨두주파수에 대한 위상으로 치환하는 단계는,

   지수감쇠코사인함수

   『

   x(t)=Ae^-αtcos(ω₁t+φ₁)

   』의

   이산푸리에 변환에서 첨두주파수 ω₁에 대응하는 위상 φ₁은

   『

   φ₁=∠X(ω₁)

   』과 같이

   이산치 푸리에 스펙트럼의 위상(∠X(ω₁))을 시계열 데이터에 포함된 첨두주파수에 대한 위상(φ₁)으로 치환하는 스텝임을 특징으로 하는 이산푸리에 변환에 의한 시계열 데이터의 위상 추정 방법.

   단, 지수감쇠코사인함수에서, ω₁은 주파수, A는 크기, α는 감쇠정수, φ₁은 첨두주파수 ω₁에 대응하는 위상이다.
3. 청구항 1에 있어서,

   상기 이산치 푸리에스펙트럼의 위상을 시계열 데이터에 포함된 첨두주파수에 대한 위상으로 치환하는 단계는,

   지수감쇠사인함수

   『

   x(t)=Ae^-αtsin(ω₁t+φ₂)

   』의

   이산푸리에 변환에서 첨두주파수 ω₁에 대응하는 위상 φ₂는

   『

   φ₂=∠(X(ω₁)+π/2)

   』와 같이

   이산치 푸리에 스펙트럼의 위상(∠X(ω₁)+π/2)을 시계열 데이터에 포함된 첨두주파수에 대한 위상(φ₂)으로 치환하는 스텝임을 특징으로 하는 이산푸리에 변환에 의한 시계열 데이터의 위상 추정 방법.

   단, 지수감쇠사인함수에서, ω₁은 주파수, A는 크기, α는 감쇠정수, φ₂는 첨두주파수 ω₁에 대응하는 위상이다.
4. 이산 푸리에 변환에 의한 시계열 데이터의 위상 추정방법에 있어서,

   데이터취득 장치로부터 이산치 시계열 데이터를 취득하는 스텝과,

   상기 취득된 이산치 시계열 데이터에 대해서 이산푸리에변환을 수행하여 이산치 푸리에 스펙트럼을 계산하여, 첨두치를 선택하는 스텝과,

   상기 첨두치의 좌측과 우측 푸리에스펙트럼을 추출하여 좌우측스펙트럼의 편차크기(ΔX₁₂)와 첨두치와 좌우측스펙트럼 중에서 작은 값의 편차크기(ΔX₁)를 계산하는 스텝과,

   상기 계산된 편차크기 비(ΔX₁₂/X₁)의 백분율을 계산하여, 편차크기 비가 10%이하면 정확한 위상으로 판정하는 스텝과,

   이 판정 결과를 출력하는 스텝을 포함하는 것을 특징으로 하는 이산푸리에 변환에 의한 시계열 데이터의 위상 추정방법.

Images