주파수 영역 분석에 대한 실용적 입문
이 예제에서는 기본적인 주파수 영역 신호 분석을 수행하고 해석하는 방법을 보여줍니다. 이 예제에서는 신호의 주파수 영역 표현을 사용하는 경우와 시간 영역 표현을 사용하는 경우의 이점을 비교하여 설명하고 시뮬레이션된 데이터와 실제 데이터 사용하여 신호의 기본 개념을 다룹니다. 또한 'FFT의 크기와 위상은 무엇을 의미합니까?', '내 신호가 주기적입니까?', '전력은 어떻게 측정합니까?', '이 대역에 신호가 하나 또는 둘 이상이 있습니까?' 등의 기본적인 질문에 대한 답을 제시합니다.
주파수 영역 분석은 신호 처리 응용에서 가장 중요한 툴입니다. 주파수 영역 분석은 통신, 지질학, 원격탐사, 영상 처리와 같은 영역에서 널리 사용됩니다. 시간 영역 분석은 신호가 시간이 지남에 따라 어떻게 변하는지를 보여주는 반면, 주파수 영역 분석은 신호의 에너지가 주파수 범위에서 어떻게 분포되는지를 보여줍니다. 주파수 영역 표현에는 원래 시간 신호를 복원하기 위해 각각의 주파수 성분에 다 적용되어야 하는 위상 변위에 대한 정보도 포함됩니다.
신호는 변환이라고 하는 수학 연산자 쌍을 사용하여 시간 영역과 주파수 영역 사이에서 변환될 수 있습니다. 이에 대한 예로 푸리에 변환을 들 수 있으며, 이 변환은 함수를 여러(무한대일 수 있음) 사인파 주파수 성분의 합으로 분해합니다. 주파수 성분의 스펙트럼은 신호의 주파수 영역 표현입니다. 푸리에 역변환은 주파수 영역 함수를 시간 함수로 다시 변환합니다. MATLAB®의 fft
함수와 ifft
함수를 사용하면 신호의 이산 푸리에 변환(DFT)과 이 변환의 역을 각각 계산할 수 있습니다.
FFT의 크기와 위상 정보
신호의 주파수 영역 표현에는 각 주파수에서의 신호 크기 및 위상에 대한 정보가 들어 있습니다. 이것이 FFT 계산의 출력값이 복소수인 이유입니다. 복소수 는 가 성립하는 실수부 과 허수부 를 가집니다. 의 크기는 으로 계산되고, 의 위상은 로 계산됩니다. MATLAB 함수 abs
와 angle
을 사용하여 임의 복소수의 크기 및 위상을 구할 수 있습니다.
오디오 예제를 사용하여 신호의 크기와 위상이 어떠한 정보를 포함하는지를 파악할 수 있습니다. 이를 수행하려면 15초간의 어쿠스틱 기타 음악이 든 오디오 파일을 불러오십시오. 오디오 신호의 샘플 레이트는 44.1kHz입니다.
Fs = 44100;
y = audioread("guitartune.wav");
fft
를 사용하여 신호의 주파수 성분을 확인합니다.
NFFT = length(y); Y = fft(y,NFFT); F = ((0:1/NFFT:1-1/NFFT)*Fs).';
FFT의 출력값은 신호의 주파수 성분에 대한 정보를 포함하는 복소수 벡터입니다. 크기는 다른 성분에 비해 상대적인 주파수 성분의 강도를 표시합니다. 위상은 모든 주파수 성분이 시간상 어떻게 정렬되는지를 표시합니다.
신호에 대한 주파수 스펙트럼의 크기와 위상 성분을 플로팅합니다. 크기는 로그 스케일(dB)로 간편하게 플로팅됩니다. 주파수의 연속 함수를 볼 수 있도록 위상은 unwrap
함수를 사용하여 펼쳐집니다.
magnitudeY = abs(Y); % FFT magnitude phaseY = unwrap(angle(Y)); % FFT phase helperFrequencyAnalysisPlot1(F,magnitudeY,phaseY,NFFT)
푸리에 역변환을 주파수 영역 벡터 Y
에 적용하여 시간 신호를 복원할 수 있습니다. "symmetric"
플래그는 ifft
에 실수 값을 갖는 시간 신호를 처리하고 있음을 알립니다. 그러면 이 함수가 계산에서 수치상의 부정확성으로 인해 역변환에 나타나는 작은 허수 성분을 0으로 처리합니다. 원래 시간 신호 y
와 복원된 신호 y1
이 거의 같다는 것을 알 수 있습니다(해당 차이의 노름이 대략 1e-14
임). 이러한 두 신호 간의 매우 작은 차이는 위에서 언급한 수치상의 정확도로 인한 것이기도 합니다. 푸리에 역변환 신호 y1
을 재생하고 들어봅니다.
y1 = ifft(Y,NFFT,"symmetric");
norm(y-y1)
ans = 3.7851e-14
hplayer = audioplayer(y1,Fs); play(hplayer);
신호의 크기 응답 변경으로 생기는 결과를 보려면, 크기를 0으로 놓아 FFT 출력값에서 1kH보다 높은 주파수 성분을 직접 제거한 후 오디오 파일의 소리에 미치는 결과를 들어봅니다. 신호의 고주파 성분을 제거하는 것을 저역통과 필터링이라고 합니다.
Ylp = Y; Ylp(F>=1000 & F<=Fs-1000) = 0; helperFrequencyAnalysisPlot1(F,abs(Ylp),unwrap(angle(Ylp)), ... NFFT,"Frequency components above 1 kHz have been zeroed")
ifft
를 사용하여 필터링된 신호를 시간 영역 신호로 다시 돌립니다.
ylp = ifft(Ylp,"symmetric");
신호를 재생합니다. 여전히 멜로디를 들을 수는 있지만 귀를 막고 듣는 것처럼 들립니다(이를 수행하면 고주파음이 필터링됨). 기타의 현을 퉁길 때 기타가 400kHz와 1kHz 사이의 음을 생성하더라도 현은 기본 주파수의 배수로 진동합니다. 고조파라고 하는 이러한 더 높은 주파수 성분은 기타에 기타 특유의 음색을 주게 됩니다. 이러한 성분을 제거하면 소리가 "불투명"하게 됩니다.
hplayer = audioplayer(ylp, Fs); play(hplayer);
신호의 위상에는 노래의 음이 시간상 언제 나타나는지에 대한 중요한 정보가 포함되어 있습니다. 오디오 신호에서 위상의 중요성을 설명하기 위해 각 주파수 성분의 크기를 사용하여 위상 정보를 완전하게 제거해 보겠습니다. 참고로, 이를 수행하면 크기 응답이 변경되지 않은 상태로 유지됩니다.
% Take the magnitude of each FFT component of the signal Yzp = abs(Y); helperFrequencyAnalysisPlot1(F,abs(Yzp), ... unwrap(angle(Yzp)),NFFT,[],"Phase has been set to zero")
시간 영역으로 신호를 되돌리고 오디오를 재생합니다. 원래의 소리를 전혀 인식할 수 없습니다. 크기 응답이 동일하며, 이번에는 주파수 성분이 제거되지 않았지만, 음의 순서가 완벽하게 사라졌습니다. 이제 신호는 시간 0에 모두 정렬된 정현파의 무리로 구성되었습니다. 일반적으로, 필터링으로 인한 위상 왜곡은 인식할 수 없는 정도로까지 신호를 손상시킬 수 있습니다.
yzp = ifft(Yzp,"symmetric");
hplayer = audioplayer(yzp, Fs);
play(hplayer);
신호의 주기성 찾기
신호의 주파수 영역 표현을 통해, 시간 영역에서 신호를 확인할 때 쉽게 볼 수 없거나 전혀 볼 수 없는 여러 가지 신호 특성을 확인할 수 있습니다. 예를 들어, 신호의 주기적 동작을 찾으려는 경우 주파수 영역 분석이 유용합니다.
사무실 건물에서 온도의 주기적 동작 분석
겨울 동안 사무실 건물에서 측정한 일련의 온도 값이 있다고 가정하겠습니다. 온도 값은 약 16.5주 동안 30분 간격으로 측정되었습니다. 주 단위로 스케일링된 시간 축을 갖는 시간 영역 데이터를 살펴보겠습니다. 이 데이터에 주기적 동작이 있을 수 있을까요?
load officetemp.mat Fs = 1/(60*30); % Sample rate is 1 sample every 30 minutes t = (0:length(temp)-1)/Fs; helperFrequencyAnalysisPlot2(t/(60*60*24*7),temp, ... "Time (weeks)","Temperature (\circF)")
시간 영역 신호를 봐서는 사무실 온도에 주기적 동작이 있는지를 파악하는 것이 거의 불가능합니다. 그러나, 주파수 영역 표현을 보면 온도의 주기적 동작이 분명하게 드러납니다.
신호의 주파수 영역 표현을 구합니다. 주당 주기 수 단위로 스케일링된 주파수 축을 갖는 FFT 출력값의 크기를 플로팅하면 다른 주파수 성분보다 명확하게 큰 2개의 스펙트럼 선이 있는 것을 확인할 수 있습니다. 한 스펙트럼 선은 1주기/주에 있고, 다른 스펙트럼 선은 7주기/주에 있습니다. 이는 온도가 조절되는 건물에서 일주일 동안 측정된 데이터를 가져온 것이란 것을 생각해보면 이해할 수 있습니다. 첫 번째 스펙트럼 선은 건물 온도가 주말 동안 낮았다가 주중에 높아지는, 주별 주기를 따른다는 것을 나타냅니다. 두 번째 스펙트럼 선은 야간에는 온도가 낮았다가 주간에는 온도가 높아지는 일일 주기도 있음을 나타냅니다.
NFFT = length(temp); % Number of FFT points F = (0 : 1/NFFT : 1/2-1/NFFT)*Fs; % Frequency vector TEMP = fft(temp,NFFT); TEMP(1) = 0; % remove the DC component for better visualization helperFrequencyAnalysisPlot2(F*60*60*24*7,abs(TEMP(1:NFFT/2)), ... "Frequency (cycles/week)","Magnitude",[],[],[0 10])
전력 측정하기
periodogram
함수는 신호의 FFT를 계산하고 출력값을 정규화하여 파워 스펙트럼 밀도(PSD)나 파워 스펙트럼을 구합니다. 이 파워 스펙트럼을 통해 전력을 측정할 수 있습니다. PSD는 시간 신호의 전력이 주파수에 대해 어떻게 분포되는지를 설명하며, 단위는 와트/Hz입니다. PSD의 각각의 점들이 정의된 주파수 구간(예: PSD의 분해능 대역폭)에 대해 그 점들을 적분하여 파워 스펙트럼을 계산합니다. 파워 스펙트럼의 단위는 와트입니다. 구간에 대해 적분할 필요 없이 파워 스펙트럼에서 직접 전력 값을 읽을 수 있습니다. 참고로, PSD와 파워 스펙트럼은 실수이므로, 위상 정보가 포함되어 있지 않습니다.
비선형 전력 증폭기의 출력값에서 고조파 측정
형태의 3차 왜곡을 가진 전력 증폭기 출력에서 측정한 데이터를 불러옵니다. 여기서 는 출력 전압이고, 는 입력 전압입니다. 데이터는 3.6kHz의 샘플 레이트로 캡처되었습니다. 입력값 는 단위 진폭을 갖는 60Hz 정현파로 구성되어 있습니다. 비선형 왜곡의 특성으로 인해 증폭기 출력 신호에는 DC 성분, 60Hz 성분, 그리고 120Hz와 180Hz에서 2차 고조파 및 3차 고조파가 포함된다고 예상할 수 있습니다.
증폭기 출력값의 3600개 샘플을 불러오고, 파워 스펙트럼을 계산하여 그 결과를 로그 스케일(데시벨-와트, 즉 dBW)로 플로팅합니다.
load ampoutput1.mat Fs = 3600; NFFT = length(y); % Power spectrum is computed when you pass a "power" flag input [P,F] = periodogram(y,[],NFFT,Fs,"power"); helperFrequencyAnalysisPlot2(F,10*log10(P),"Frequency (Hz)", ... "Power spectrum (dBW)",[],[],[-0.5 200])
파워 스펙트럼의 플롯은 예상되었던 4개 피크 중 3개를 DC, 60Hz, 120Hz에서 보여줍니다. 또한, 신호의 잡음으로 인해 발생하는 스퓨리어스 피크도 몇 개 더 보여줍니다. 180Hz 고조파는 잡음에 완전히 묻혀 있음을 알 수 있습니다.
표시되는 예상 피크의 전력을 측정합니다.
PdBW = 10*log10(P); power_at_DC_dBW = PdBW(F==0)
power_at_DC_dBW = -7.8873
[peakPowers_dBW, peakFreqIdx] = findpeaks(PdBW,MinPeakHeight=-11); peakFreqs_Hz = F(peakFreqIdx)
peakFreqs_Hz = 2×1
60
120
peakPowers_dBW
peakPowers_dBW = 2×1
-0.3175
-10.2547
잡음이 있는 신호의 전력 측정값 향상
위의 플롯에서 볼 수 있듯이, 주기도는 관심 신호와 상관없는 여러 주파수 피크를 보여줍니다. 스펙트럼에 잡음이 매우 많아 보입니다. 이렇게 보이는 이유는 잡음이 있는 신호에 대해 하나의 짧은 구현만 분석했기 때문입니다. 실험을 여러 번 반복하고 평균을 구하면 스퓨리어스 스펙트럼 피크가 제거되고 더욱 정확한 전력 측정값이 생성됩니다. pwelch
함수를 사용하여 이러한 평균화 작업을 수행할 수 있습니다. 이 함수는 큰 데이터 벡터를 받아서 지정된 길이의 더 작은 세그먼트로 분할하고, 세그먼트 개수와 동일한 개수의 주기도를 계산하고, 이에 대한 평균을 구합니다. 사용 가능한 세그먼트 수가 늘어나면 pwelch
함수가 예상 값에 더 가까운 전력 값을 사용하여 더 평활화된 파워 스펙트럼(분산이 더 적음)을 생성합니다.
증폭기 출력값의 500e3
개 점으로 구성된 더 큰 규모의 관측값을 불러옵니다. floor(500e3/3600)
= 138 FFT의 평균으로 파워 스펙트럼을 구할 수 있도록, FFT를 수행하는 데 사용된 점의 개수를 3600
개로 유지합니다.
load ampoutput2.mat SegmentLength = NFFT; % Power spectrum is computed when you pass a "power" flag input [P,F] = pwelch(y,ones(SegmentLength,1),0,NFFT,Fs,"power"); helperFrequencyAnalysisPlot2(F,10*log10(P),"Frequency (Hz)", ... "Power spectrum (dBW)",[],[],[-0.5 200])
플롯에서 볼 수 있듯이, pwelch
는 잡음으로 인해 발생하는 모든 스퓨리어스 주파수 피크를 효과적으로 제거합니다. 잡음에 묻혀 있었던 180Hz의 스펙트럼 성분이 이제 표시됩니다. 평균화를 수행하면 스펙트럼에서 분산이 제거되어 더 정확한 전력 측정값을 얻을 수 있습니다.
총 평균 전력과 특정 주파수 대역에 대한 총 전력 측정
시간 영역 신호의 총 평균 전력을 측정하는 것은 흔하고 손쉬운 작업입니다. 증폭기 출력 신호 y
에 대해 시간 영역에서 총 평균 전력은 다음과 같이 계산됩니다.
pwr = sum(y.^2)/length(y) % in watts
pwr = 8.1697
주파수 영역에서는 총 평균 전력이 신호 내 모든 주파수 성분에 대한 전력의 합으로 계산됩니다. pwr1
의 값은 신호의 파워 스펙트럼에서 사용할 수 있는 모든 주파수 성분의 합으로 구성됩니다. 이 값은 시간 영역 신호를 사용하여 위에서 계산된 pwr
값과 일치합니다.
pwr1 = sum(P) % in watts
pwr1 = 8.1698
하지만, 주파수의 한 대역에서 사용할 수 있는 총 전력을 측정하려면 어떻게 해야 할까요? bandpower
함수를 사용하여 원하는 주파수 대역에 대해 전력을 계산할 수 있습니다. 시간 영역 신호를 이 함수에 대한 입력값으로 직접 전달하여 지정된 대역에 대한 전력을 구할 수 있습니다. 이 경우, 이 함수는 주기도를 사용해 파워 스펙트럼을 추정합니다.
50Hz~70Hz 대역에서 전력을 계산합니다. 결과값에는 60Hz 전력과 관심 대역에 대한 잡음 전력이 포함됩니다.
pwr_band = bandpower(y,Fs,[50 70]);
pwr_band_dBW = 10*log10(pwr_band) % dBW
pwr_band_dBW = 0.0341
특정 대역에서 전력을 측정하는 데 사용된 파워 스펙트럼의 계산을 제어하려는 경우 PSD 벡터를 bandpower
함수에 전달할 수 있습니다. 예를 들어, 앞에서 한 것처럼 pwelch
함수를 사용하여 PSD를 계산하고 잡음 효과의 평균화를 수행할 수 있습니다.
% Power spectral density is computed when you specify the "psd" option [PSD,F] = pwelch(y,ones(SegmentLength,1),0,NFFT,Fs,"psd"); pwr_band1 = bandpower(PSD,F,[50 70],"psd"); pwr_band_dBW1 = 10*log10(pwr_band1) % dBW
pwr_band_dBW1 = 0.0798
스펙트럼 성분 찾기
신호는 하나 이상의 주파수 성분으로 구성될 수 있습니다. 모든 스펙트럼 성분을 확인할 수 있는 능력은 분석의 주파수 분해능에 따라 달라집니다. 주파수 분해능이나 파워 스펙트럼의 분해능 대역폭은 R = Fs/N으로 정의되며, 여기서 N은 신호 관측값의 길이입니다. 주파수 분해능보다 큰 주파수로 분리된 스펙트럼 성분만 분해됩니다.
건물의 지진 진동 제어 시스템 분석
능동질량감쇠기(AMD) 제어 시스템은 지진이 발생했을 때 건물의 진동을 줄이는 데 사용됩니다. 능동질량감쇠기는 건물의 최상층에 설치되며, 건물 각 층의 변위 및 가속도 측정값에 따라 제어 시스템이 신호를 감쇠기에 전송해 지반 교란이 감쇠되도록 질량을 이동시킵니다. 지진 조건하에서 3층짜리 테스트 구조물의 1층에서 가속도 측정값이 기록되었습니다. 가속도 값은 능동질량감쇠기 제어 시스템이 없는 상태(개루프 조건)와 능동질량감쇠기 제어 시스템이 있는 상태(폐루프 조건)에서 측정되었습니다.
가속도 데이터를 불러오고 1층의 가속도에 대해 파워 스펙트럼을 계산합니다. 데이터 벡터의 길이는 10e3
이고 샘플 레이트는 1kHz입니다. 64개 데이터 점의 길이를 가지는 세그먼트에 pwelch
를 사용하여 FFT 평균값(floor(10e3/64)
= 156)과 분해능 대역폭(Fs/64
= 15.625Hz)을 구합니다. 위에서 설명한 바와 같이, 평균화를 수행하면 잡음 효과가 줄어들고 더 정확한 전력 측정값을 얻을 수 있습니다. 512개 FFT 점을 사용합니다. NFFT > N을 사용하면 주파수 점이 보간되게 되어 더 자세한 스펙트럼 플롯이 생성됩니다. 이는 시간 신호의 끝부분에 NFFT-N개의 0을 추가하고 0으로 채워진 벡터의 NFFT개 점을 갖는 FFT를 구하는 방식으로 수행됩니다.
개루프와 폐루프의 가속도 파워 스펙트럼을 보면 제어 시스템이 활성 상태일 때 가속도 파워 스펙트럼이 4dB에서 11dB까지 줄어드는 것을 알 수 있습니다. 최대 감쇠는 약 23.44kHz에서 발생합니다. 11dB 감소는 진동 전력이 12.6배 감소됨을 의미합니다. 총 전력은 0.1670와트에서 0.059와트로 2.83배 감소됩니다.
load quakevibration.mat Fs = 1e3; % sample rate NFFT = 512; % number of FFT points segmentLength = 64; % segment length % open loop acceleration power spectrum [P1_OL,F] = pwelch(gfloor1OL,ones(segmentLength,1),0,NFFT,Fs,"power"); % closed loop acceleration power spectrum P1_CL = pwelch(gfloor1CL,ones(segmentLength,1),0,NFFT,Fs,"power"); helperFrequencyAnalysisPlot2(F,10*log10([(P1_OL) (P1_CL)]),... "Frequency (Hz)","Acceleration Power Spectrum (dB)",... "Resolution bandwidth = 15.625 Hz",["Open" "Closed"]+" loop",[0 100])
진동 데이터의 분석을 통해 진동이 주기적 동작을 가진다는 것을 알 수 있습니다. 그러면 왜 위에 보이는 스펙트럼 플롯에는 일반적으로 주기적 동작에서 볼 수 있는 예리한 스펙트럼 선이 포함되어 있지 않을까요? 64개 점 세그먼트 길이로 구한 분해능에서는 이러한 선을 분해할 수 없기 때문에 누락되었던 건 아닐까요? 이전에 분해할 수 없었던 스펙트럼 선이 있는지 확인하기 위해 주파수 분해능을 늘려보십시오. pwelch
함수에 사용된 데이터 세그먼트의 길이를 512개 점으로 늘려서 이 작업을 수행합니다. 그러면 새 분해능(Fs/512
= 1.9531Hz)이 생성됩니다. 이 경우, FFT 평균 개수가 floor(10e3/512)
= 19로 줄어듭니다. pwelch
를 사용하는 경우 평균 개수와 주파수 분해능 간에는 명백한 상호 절충 관계가 있습니다. FFT 점의 개수는 512로 그대로 유지합니다.
NFFT = 512; % number of FFT points segmentLength = 512; % segment length [P1_OL,F] = pwelch(gfloor1OL,ones(segmentLength,1),0,NFFT,Fs,"power"); P1_CL = pwelch(gfloor1CL,ones(segmentLength,1),0,NFFT,Fs,"power"); helperFrequencyAnalysisPlot2(F,10*log10([(P1_OL) (P1_CL)]), ... "Frequency in Hz","Acceleration Power Spectrum in dB", ... "Resolution bandwidth = 1.95 Hz",{"Open loop", "Closed loop"},[0 100])
주파수 분해능을 늘린 결과 개루프 스펙트럼에서는 3개 피크가 관찰되고 폐루프 스펙트럼에서는 2개 피크가 관찰됨을 알 수 있습니다. 이러한 피크는 이전에는 분해할 수 없었습니다. 개루프 스펙트럼에서 피크 간 간격은 약 11Hz로, 길이가 64인 세그먼트에서 구한 주파수 분해능보다는 작고, 길이가 512인 세그먼트에서 구한 분해능보다는 큽니다. 이제 진동의 주기적 동작을 볼 수 있습니다. 주진동 주파수는 5.86Hz이며, 주파수 피크가 균일한 간격으로 있는 것은 서로 고조파 관계임을 나타냅니다. 제어 시스템이 진동의 전체 전력을 줄인다는 사실은 이미 확인되었지만 분해능이 더 높은 스펙트럼에서 보면 제어 시스템이 17.58Hz에서 고조파 성분을 상쇄하는 역할도 한다는 것을 알 수 있습니다. 따라서 제어 시스템은 진동을 줄일 뿐만 아니라 정현파에 더 가깝게 합니다.
주파수 분해능은 FFT 점 개수가 아니라 신호 점 개수로 결정된다는 점을 유의해야 합니다. FFT 점 개수를 늘리면 주파수 데이터가 보간되어 스펙트럼에 대한 더 세부적인 정보를 제공하지만, 분해능은 향상시키지 않습니다.
결론
이 예제에서는 fft
,
ifft
, periodogram
, pwelch
, bandpower
함수를 사용하여 신호의 주파수 영역 분석을 수행하는 방법을 배웠습니다. FFT의 복소수적 특성을 이해하고 주파수 스펙트럼의 크기와 위상에 포함된 정보를 살펴봤습니다. 신호의 주기성을 분석할 때 주파수 영역 데이터를 사용함으로써 얻을 수 있는 이점을 확인했습니다. 잡음이 있는 신호에 대해 총 전력을 계산하거나 특정 주파수 대역에서 전력을 계산하는 방법을 배웠습니다. 스펙트럼의 주파수 분해능을 늘리면 간격이 좁은 주파수 성분을 관찰할 수 있다는 것을 이해했고 주파수 분해능과 스펙트럼 평균화 간의 상호 절충에 대해서도 배웠습니다.
추가 참고 자료
주파수 영역 분석에 대한 자세한 내용은 스펙트럼 분석 항목을 참조하십시오.
참고 문헌: J.G. Proakis and D. G. Manolakis, "Digital Signal Processing. Principles, Algorithms, and Applications", Prentice Hall, 1996.
부록
이 예제에서는 다음 헬퍼 함수가 사용되었습니다.
참고 항목
fft
| periodogram
| pspectrum
| pwelch