Main Content

불균일하게 샘플링된 신호의 스펙트럼 분석

이 예제에서는 불균일하게 샘플링된 신호에 대해 스펙트럼 분석을 수행하는 방법을 보여줍니다. 신호가 불균일하게 샘플링되었는지 여부를 확인하고, 불균일하게 샘플링되지 않았다면 신호의 스펙트럼 또는 파워 스펙트럼 밀도를 계산하는 방법을 보여줍니다.

이 예제에서는 불균일하게 샘플링된 신호의 스펙트럼을 계산할 수 있는 Lomb-Scargle 주기도를 소개합니다.

불균일하게 샘플링된 신호

불균일하게 샘플링된 신호는 자동차 산업에서, 통신에서, 그리고 의료와 천문학을 비롯한 다양한 분야에서 찾을 수 있습니다. 불균일한 샘플링은 불완전한 센서, 불일치된 클록 또는 이벤트로 트리거되는 현상으로 인해 발생할 수 있습니다.

스펙트럼 성분의 계산과 연구는 신호 분석의 중요한 부분입니다. 주기도, Welch 방법과 같은 일반적인 스펙트럼 분석 기법에는 균일하게 샘플링된 입력 신호가 필요합니다. 샘플링이 균일하지 않은 경우, 신호를 리샘플링하거나 균일한 샘플 그리드로 보간할 수 있습니다. 그러나 이로 인해 스펙트럼에 원치 않는 아티팩트가 추가되어 분석 오류가 발생할 수 있습니다.

이보다 나은 대안은 불균일한 샘플을 직접 사용하여 리샘플링이나 보간이 필요하지 않은 Lomb-Scargle 방법을 사용하는 것입니다. 이 알고리즘은 plomb 함수에 구현되어 있습니다.

누락된 데이터가 있는 신호에 대한 스펙트럼 분석

마이크로컨트롤러가 방의 온도를 기록하고 측정값이 저장되는 클라우드 기반 서버로 15분에 한 번씩 측정값을 전송하는 온도 모니터링 시스템이 있다고 가정하겠습니다. 이 시스템에는 인터넷 연결에 결함이 있어 마이크로컨트롤러가 전송하는 측정값의 일부를 클라우드 기반 시스템이 수신할 수 없습니다. 또한, 측정 기간 동안 마이크로컨트롤러의 배터리가 최소 한 번 방전되어 샘플링에 큰 공백이 발생했습니다.

온도 측정값과 이에 대응하는 타임스탬프를 불러오겠습니다.

load('nonuniformdata.mat','roomtemp','t1')

figure
plot(t1/(60*60*24*7),roomtemp,'LineWidth',1.2)

grid
xlabel('Time (weeks)')
ylabel('Temperature (\circF)')

Figure contains an axes object. The axes object with xlabel Time (weeks), ylabel Temperature ( degree F) contains an object of type line.

신호가 균일하게 샘플링되었는지 여부를 확인하는 간단한 방법은 연속된 샘플 시간 사이의 간격을 히스토그램으로 살펴보는 것입니다.

샘플링 간격(시간 차)의 히스토그램을 분 단위로 플로팅합니다. 샘플이 존재하는 점만 포함하여 플로팅하십시오.

tAtPoints = t1(~isnan(roomtemp))/60;
TimeIntervalDiff = diff(tAtPoints);

figure
hist(TimeIntervalDiff,0:100)
grid
xlabel('Sampling intervals (minutes)')
ylabel('Occurrences')
xlim([10 100])

Figure contains an axes object. The axes object with xlabel Sampling intervals (minutes), ylabel Occurrences contains an object of type patch. This object represents TimeIntervalDiff.

측정값의 대부분이 예상대로 약 15분의 간격을 갖습니다. 그러나 약 30분과 45분의 간격을 갖는 측정값도 꽤 있는데, 이는 샘플이 1회 또는 2회 연속하여 누락된 경우에 해당합니다. 이로 인해 신호가 불균일하게 샘플링되었습니다. 또한, 히스토그램에서 높은 빈도를 나타내는 막대 주변에서 얼마간의 지터를 볼 수 있습니다. 이것은 TCP/IP 대기 시간과 관련이 있을 수 있습니다.

Lomb-Scargle 방법을 사용하여 신호의 스펙트럼 성분을 계산하고 시각화합니다. 스펙트럼을 보다 잘 시각화할 수 있도록 주당 약 13번의 사이클에 해당하는 0.02mHz까지의 주파수를 살펴봅니다.

[Plomb,flomb] = plomb(roomtemp,t1,2e-5,'power');

figure
plot(flomb*60*60*24*7,Plomb)
grid
xlabel('Frequency (cycles/week)')
ylabel('Power (dBW)')

Figure contains an axes object. The axes object with xlabel Frequency (cycles/week), ylabel Power (dBW) contains an object of type line.

이 스펙트럼은 주당 7번의 사이클과 주당 1번의 사이클에서 강한 주기성을 보입니다. 이는 온도가 조절되는 건물에서 한 주를 기준으로 하여 측정된 데이터를 가져온 것이란 것을 생각해보면 이해할 수 있습니다. 주당 1번의 사이클에서 피크를 보이는 스펙트럼 선은 건물 온도가 주말 동안 낮았다가 주중에 높아지는 주별 주기를 따른다는 것을 나타냅니다. 주당 7번의 사이클에 해당하는 스펙트럼 선은 야간에는 온도가 낮았다가 주간에는 온도가 높아지는 일별 주기도 있음을 나타냅니다.

간격이 균일하지 않는 샘플을 갖는 신호에 대한 스펙트럼 분석

심박 간 시간의 생리적 변동을 나타내는 심박 변이도(HRV) 신호는 대체로 불균일하게 샘플링됩니다. 사람의 심박수는 일정하지 않기 때문입니다. HRV 신호는 심전도(ECG) 측정값에서 도출됩니다.

HRV 신호의 샘플 점은 심전도의 R-피크에서 찾을 수 있습니다. 각 점의 진폭은 연속된 R-피크 사이의 시간 차의 역으로 계산되며 두 번째 R-피크 시점에 배치됩니다.

% Load the signal, the timestamps, and the sample rate
load('nonuniformdata.mat','ecgsig','t2','Fs')

% Find the ECG peaks
[pks,locs] = findpeaks(ecgsig,Fs, ...
    'MinPeakProminence',0.3,'MinPeakHeight',0.2);

% Determine the RR intervals
RLocsInterval = diff(locs);

% Derive the HRV signal
tHRV = locs(2:end);
HRV = 1./RLocsInterval;

% Plot the signals
figure
a1 = subplot(2,1,1); 
plot(t2,ecgsig,'b',locs,pks,'*r')
grid
a2 = subplot(2,1,2);
plot(tHRV,HRV)
grid
xlabel(a2,'Time(s)')
ylabel(a1,'ECG (mV)')
ylabel(a2,'HRV (Hz)')

Figure contains 2 axes objects. Axes object 1 with ylabel ECG (mV) contains 2 objects of type line. One or more of the lines displays its values using only markers Axes object 2 with xlabel Time(s), ylabel HRV (Hz) contains an object of type line.

R-피크 간의 변화하는 간격으로 인해 HRV 데이터의 샘플 시간 불균일성이 발생합니다. 신호의 피크 위치를 살펴보고 피크 간 간격의 히스토그램을 초 단위로 플로팅합니다.

figure
hist(RLocsInterval)

grid
xlabel('Sampling interval (s)')
ylabel('RR distribution')

Figure contains an axes object. The axes object with xlabel Sampling interval (s), ylabel RR distribution contains an object of type patch. This object represents RLocsInterval.

HRV 스펙트럼의 일반적인 관심 주파수 대역은 다음과 같습니다.

  • 초저주파수(VLF): 3.3 ~ 40mHz

  • 저주파수(LF): 40 ~ 150mHz

  • 고주파수(HF): 150 ~ 400mHz

이러한 대역은 HRV에 기여하는 상이한 생물학적 조절 메커니즘의 주파수 대역을 개략적으로 제한합니다. 이러한 대역에서 발생하는 변동은 생물학적으로 유의미합니다.

plomb를 사용하여 HRV 신호의 스펙트럼을 계산합니다.

figure
plomb(HRV,tHRV,'Pd',[0.95, 0.5])

Figure contains 2 axes objects. Axes object 1 is empty. Axes object 2 with title Lomb-Scargle Power Spectral Density Estimate, xlabel Frequency (mHz), ylabel Power/frequency (dB/Hz) contains 3 objects of type line.

파선은 95% 및 50% 검출 확률을 나타냅니다. 이러한 임계값은 피크의 통계적 유의성을 측정합니다. 스펙트럼은 위에 나열된 세 개의 관심 대역 모두에서 피크를 보입니다. 그러나 VLF 범위에서 23.2mHz에 위치한 피크만 검출 확률 95%를 보이고, 나머지 피크는 검출 확률이 50%보다 낮습니다. 40mHz 아래에 있는 피크는 체온 조절 시스템과 호르몬 요인과 같은 장기적 조절 메커니즘의 결과로 간주됩니다.

참고 항목