실습 C. 표본화, 창, 시간-주파수 분석

표본화율 f_s는 보이는 주파수 범위, 표본 수 N은 DFT 칸 간격, 창은 누설을 정한다.

이 실습의 자리: 시뮬레이션 — 실습 B의 DFT를 바탕으로 표본화와 창의 영향을 보고, 다음 실습 D에서 비트를 파형으로 보내고 되찾는다.

준비

활동은 셋이다. 접힌 주파수를 구하는 함수를 쓰고 그림과 비교한 뒤, 표본 수와 DFT 길이 손잡이, 창 선택으로 스펙트럼을 관찰한다. 모든 입력은 합성 신호이며 Pluto는 쓰지 않는다. 셀은 실습 A의 준비처럼 위에서부터 실행한다.

이론은 2장 표본화와 aliasing, 3장 표본 수와 주파수 간격, 3장 누설과 창 함수를 따른다.

그림과 표 옆의 기록에 넣기로 결과를 담고, 글 상자의 설명과 함께 도구 막대의 저장으로 저장한다.

아래 도우미는 창을 곱한 기록의 DFT를 구한다. 창의 합으로 나누므로 칸 위에 놓인 단위 진폭 복소 톤을 기준으로 크기를 읽는다. fft_length가 기록보다 길면 뒤에 영을 덧붙인다.

from sdrlab import simviews


def spectrum_db(x, sample_rate, fft_length=None, window='rect'):
    """신호 x와 표본화율(Hz)을 받아 주파수(Hz), 창 보정 크기(dB)를 돌려준다."""
    x = np.asarray(x)
    window_weights = show.window(window, x.size)
    fft_length = fft_length or x.size

    spectrum = np.fft.fftshift(np.fft.fft(x * window_weights, fft_length))
    f = np.fft.fftshift(np.fft.fftfreq(fft_length, 1 / sample_rate))
    magnitude = np.abs(spectrum) / window_weights.sum()
    magnitude_db = 20 * np.log10(np.maximum(magnitude, 1e-12))
    return f, magnitude_db

활동 1. 대역 밖의 주파수 접기

표본 사이에서 더 돈 바퀴 수를 알 수 있을까? $f$와 $f+mf_s$의 복소 톤은 정수 $m$에 대해 같은 표본을 만든다. 그래서 대표 주파수를 $-f_s/2\le f<f_s/2$에 고른다. 실수 코사인은 양·음 회전을 함께 담으므로 접힌 주파수의 양쪽에 대칭 성분이 있다. 복소 톤은 그 방향까지 구별한다.

쓸 함수 fold(f, sample_rate)

• 입력: 주파수 f (Hz, 수 하나 또는 배열), 양수인 표본화율 sample_rate (Hz)
• 출력: 같은 모양의 접힌 주파수 (Hz), 범위는 $[-f_s/2,f_s/2)$

계산 순서

1. np.asarray(f)로 입력을 배열로 바꾸고 sample_rate / 2를 더한다.
2. np.remainder(..., sample_rate)로 표본화율로 나눈 나머지를 구한다. 양수인 sample_rate의 나머지는 음수가 아니다.
3. sample_rate / 2를 빼서 중심을 원점으로 되돌린다.

def fold(f, sample_rate):
    """입력 f(Hz), sample_rate(Hz)에서 같은 표본을 만드는 접힌 주파수(Hz)를 돌려준다."""
    # TODO: 위의 계산 순서대로 fold 함수를 구현하세요.
    raise NotImplementedError('fold 함수를 구현하세요.')

아래 셀은 표본화율 1 MHz에서 입력 300, 700, 1300, 2450, −600 kHz를 비교한다. 접힌 주파수는 차례로 300, −300, 300, 450, 400 kHz이며, 복소 톤 FFT의 봉우리와 같아야 한다. 경계에서는 +500 kHz도 −500 kHz로 표시한다.

sample_rate = 1e6  # Hz
sample_count = 1000
input_frequencies = np.array([300e3, 700e3, 1300e3, 2450e3, -600e3])  # Hz

folded_frequencies = fold(input_frequencies, sample_rate)
measured_peaks = simviews.compare_fold(
    input_frequencies,
    folded_frequencies,
    sample_rate,
    sample_count,
    spectrum_db,
)

correct = (
    np.allclose(folded_frequencies, measured_peaks)
    and fold(sample_rate / 2, sample_rate) == -sample_rate / 2
)
print('→ 맞음' if correct else '→ 다시 본다')

이제 입력 주파수를 바꾼다. 먼저 복소 톤의 봉우리가 원점의 어느 쪽에 있는지 보고, 그다음 실수 코사인의 양쪽 봉우리를 본다. 표에서 식과 실제 봉우리를 함께 읽는다. 원점과 대역 경계에서는 실수 신호의 두 성분이 같은 자리에 겹친다.

def draw_alias(input_khz):
    simviews.alias(input_khz, sample_rate, sample_count, fold, spectrum_db)


show.control(
    draw_alias,
    input_khz=show.slider(-1500, 2500, step=50, value=700, unit='kHz', label='입력 주파수'),
)

활동 2. 표본을 더 모으기와 영 덧붙이기

DFT 점을 촘촘히 하면 봉우리 자체도 좁아질까? 관측 길이는 $T=N/f_s$, 원래 DFT 칸 간격은 $f_s/N$이다. 뒤에 영을 덧붙여 DFT 길이를 fft_length로 늘리면 표시점 간격은 $f_s/\mathrm{nfft}$가 되지만, 같은 유한 기록의 스펙트럼을 더 촘촘히 읽을 뿐이다. 가까운 톤을 가르는 봉우리 폭을 줄이려면 실제 관측을 늘려야 한다.

먼저 N을 고정하고 fft_length를 바꿔 표시점만 비교한다. 그다음 N을 바꾸고, 같은 기록을 촘촘히 계산한 배경 곡선의 폭을 비교한다. 입력은 123.4 Hz, 표본화율은 1000 Hz이다. N이 100에서 400으로 늘면 관측 길이는 0.1초에서 0.4초로 늘고 원래 칸 간격은 10 Hz에서 2.5 Hz로 줄어든다.

def draw_padding(N, fft_length):
    simviews.padding(N, fft_length, spectrum_db)


show.control(
    draw_padding,
    N=show.choice([100, 400], value=100, label='표본 수 N'),
    fft_length=show.choice([400, 1000, 10000], value=1000, label='DFT 길이'),
)

활동 3. 창으로 큰 톤의 누설 낮추기

큰 톤 옆의 작은 톤을 어떻게 드러낼까? 유한 구간을 자르면 큰 톤의 성분이 주위로 퍼지는 누설(leakage)이 생긴다. 창은 기록 양 끝을 완만하게 줄여 옆 봉우리를 낮추지만 중심 봉우리는 넓힌다. 칸 사이에 꼭대기가 놓이면 DFT가 꼭대기를 비껴 읽어 같은 진폭도 낮게 읽히므로 창과 칸 위치를 함께 보아야 한다.

큰 톤은 8.3 Hz, 작은 톤은 12.5 Hz이고 진폭은 각각 1과 0.01로, 작은 톤이 40 dB 낮다. 먼저 큰 톤만의 누설과 두 톤을 더한 선을 견주고, 그다음 창을 바꾼다. 사각 창에서는 큰 톤의 옆 봉우리가 작은 톤을 가린다. Hann에서는 두 톤의 겹침이 봉우리의 자리와 높이를 바꾼다. Blackman은 작은 톤을 제 높이 가까이 드러내지만, 남은 누설 때문에 봉우리 위치는 실제 주파수와 조금 다르다. 원래 창에 맞춰 크기를 정규화했으므로 창이 바뀌었다고 전체 진폭이 작아진 것으로 읽지 않는다.

def draw_leakage(window):
    simviews.leakage(window, spectrum_db)


show.control(
    draw_leakage,
    window=show.choice(['rect', 'hann', 'blackman'], value='rect', label='창'),
)