190527 DSP 노트

spectral leakage, coherent gain, scalloping loss

https://www.recordingblogs.com/wiki/coherent-gain

• spectral leakage: DFT시 디지털 샘플로 정확히 원하는 주파수의 magnitude를 나타낼 수 없어(정해진 정수개의 sample로 커버가 안되는 주기 같은 것) 인접 주파수에 spectrum이 번져 보이는 현상. window를 적용하면 spectral leakage가 완화됨. 주기 신호처럼 만들어주기 때문인 듯.
• coherent gain: window의 평균 값(전부 더한 값/윈도우 길이)
• window마다 coherent gain이 다르기 때문에, window를 적용하면 loss가 생기는 것을 방지 하기 위해 coherent gain으로 나눠주어 normalization을 수행함.

http://web.mit.edu/xiphmont/Public/windows.pdf https://www.recordingblogs.com/wiki/scalloping-loss https://ipfs.io/ipfs/QmXoypizjW3WknFiJnKLwHCnL72vedxjQkDDP1mXWo6uco/wiki/Windowfunction.html http://www.bores.com/courses/advanced/windows/10pl.htm

• spectrum에서 가장 큰 낙폭을 보이는 주파수 bin의 중간 지점의 attenuation. 가장 큰 낙폭은 주로 main lobe에서 일어나므로 main lobe에 대해 주로 말하는 것 같다. 그런데 이 bin 사이의 지점에 대한 정보는 어떻게 얻는다는 것인지 이해가 안된다… 그냥 이 개념이 이해가 안됨.

zero-padding

https://dsp-nbsphinx.readthedocs.io/en/nbsphinx-experiment/spectralanalysisdeterministicsignals/zeropadding.html#Relation-between-Discrete-Fourier-Transformations-with-and-without-Zero-Padding

• 예를 들자면 [1, 2, 3, 1] 이라는 time domain 신호를 [1, 2, 3, 1, 0, 0, 0, 0]으로 만들어 원래 N(=4) 포인트 DFT를 하던 것을 M(=8) DFT 하겠다는 것
• 이렇게 하면 주파수 상에서 기존 N 포인트로 보던 스펙트럼을 M 포인트로 볼 수 있으니 더 해상도가 높아진 것과 같은 효과가 남
• 하지만 해상도가 높아졌다고 해서 없던 정보를 더 얻은 것은 아니고 단순히 보기 좋아지게끔 interpolation 된 것. 용도에 따라 쓸모 있는 기법이 될 수 있음.

zero-padding과 DFT, DTFT의 관계

• DTFT는 time domain 상에서만 discrete이기 때문에 freq. domain에서는 연속임. 이를 컴퓨터 상에터 보려면 주파수 도메인 bin을 아주 촘촘히 만들어 연속인 것 같이 나타내면 됨.
• 이 (이산이지만 아주 촘촘한) DTFT를 사실은 DFT로부터 만들어낼 수도 있는데, DFT에 IDFT하고 다시 DTFT 하면 됨(…). 이 과정 (https://dsp-nbsphinx.readthedocs.io/en/nbsphinx-experiment/spectralanalysisdeterministicsignals/zeropadding.html#Interpolation-of-the-Discrete-Fourier-Transformation) 에서 수식에 sum이 두 번 들어가는데 periodic sinc function을 쓰면 sum을 한번만 쓴 형태로 바꾸어 나타낼 수 있음. (https://en.wikipedia.org/wiki/Dirichlet_kernel의 Variant of identity 부분을 볼 것) 처음엔 psinc로 나타낸 수식에 무슨 물리적 의미가 있나 한참 고민했는데, 정신건강에 좋게 그냥 수학적 간결함을 위한 것이라 합리화하였음. (리마인드: sinc function은 사각 윈도우의 주파수 응답. 샘플링 시 impulse 형태가 아니라 어쩔 수 없는 구형파 형태로 필터링 되기 때문에 sampling function이라고도 함 https://www.quora.com/What-is-the-difference-between-sinc-and-sampling-function )
• 위에서 말했듯이 디지털 상에서 억지로 나타낸 DTFT는 사실 원래 신호를 매우 해상도가 높은 M 포인트 DFT로 표현한 것 뿐임. 즉 수식상에서는 DFT -> IDFT -> DTFT로 만드는 과정이어도 DTFT의 표현을 디지털 상에서 M 포인트 DFT로 나타냈다는 것은 N포인트 DFT를 M 포인트 DFT로 interpolation 한 것이고, 이는 원래 time-domain 신호에 zero-padding 을 통해 M 포인트 DFT로 만든 것과 동일한 결과를 가져옴.

PSD, periodogram, welch technique

https://dsp-nbsphinx.readthedocs.io/en/nbsphinx-experiment/spectralestimationrandom_signals/periodogram.html

• PSD는 신호 고유의 값이라 알 수 있는 값이 아니라 추정할 수 있는 값임. 스펙트럼을 제곱하여 얻는 periodogram은 우리가 'PSD'로 알고 있는 것으로, PSD estimator라 보면 됨.
• periodogram을 주파수에 대해 평균내면 '주파수에 대해 평균 낸' PSD와 매우 비슷한 값이 됨.

https://dsp-nbsphinx.readthedocs.io/en/nbsphinx-experiment/spectralestimationrandomsignals/welchmethod.html

• periodogram으로 PSD를 보면 너무 들쑥날쑥하므로 STFT처럼 windowing을 하여 윈도우마다 periodogram을 계산한 후 window들에 대해 평균을 냄 (결과는 주파수마다 나오는 형태가 됨). 이것이 welch technique(method) 인데 window가 꼭 overlap이 있을 필요는 없음(overlap이 하나도 없이 window끼리 붙어 있는 형태로도 함)
• 이렇게 하면 주파수 별 PSD를 잘 반영함

causal filter, convolution, correlation

• f(t + tau)와 f(t - tau)를 이해하자. plot 상에 수평이동 되는 방향이 마치 반대로 움직이는 것 같아 헷갈리기 쉽다.
• f(t + tau): t에서 tau만큼 미래에 있는 값. f(t - tau)는 그 반대…
• 왜 수평이동은 반대 방향으로 되는 것처럼 보이는가? f(t + tau)의 plot은 '미래(t + tau)'에 있는 값을 현재 t 에 나타낸 것이기 때문에 값이 과거로 움직인 것처럼 보이는 것. 반대의 경우는 '과거'의 값을 '현재'로 가지고 오기 때문에 미래 방향으로 움직인 것처럼 보이는 것…

https://en.wikipedia.org/wiki/Convolution

• 예전부터 convolution의 수식이 너무 이해가 안됐다. 이해라기 보단 그냥 저런 공식을 쓰는 연산이다라고 생각해왔음.
• f*g(t) = intergral(f(tau)g(t - tau), dtau) 인데 여기서 f를 필터, g를 신호라고 생각하자. 이 경우 f(tau)는 필터의 tau번째 미래 tap이 되는 것이고, g(t - tau)는 t 시점으로부터 tau만큼 과거에 있는 신호값이다. 즉 f(tau)g(t - tau)란 t보다 tau 이전의 신호를 미리 f에 집어 넣었고, 그게 t 시점에 f(tau)를 만나 곱해진 값이 되는 것이다. tau 만큼 과거에 밀어넣은 값이므로, tau만큼 미래의 tap과 만나는 것은 너무나 당연하다. 이렇게 t 시점에서 tau이전부터 하나씩 넣어온 신호들을 전부 t 시점에 모아서 더하는 것이 convolution이다.

https://en.wikipedia.org/wiki/Causal_filter

• causal filter의 경우 미래의 값에 대해 필터링을 할 순 없으므로 현재시점부터 과거로만 가게 tau를 배치한다.

https://en.wikipedia.org/wiki/Cross-correlation

• cross-correlation의 경우 f, g 둘 다를 '신호'라 간주하는 것이 이해하기 쉽다. 애초 계산하는 이유도 두 신호의 similarity를 측정하기 위함.
• cross-correlation의 수식의 f, g 중 왜 한쪽은 conjugate인가? complex space에 있는 각 신호의 유사도는 스칼라여야 하기 때문에, complex 값끼리 inner product를 해야 하기 때문에 한쪽을 conjugate로 만들어 곱하고 있는 것
• cross-correlation 수식은 f(t)g(t + tau)와 f(t - tau)g(t)가 있는데, 내 눈엔 f(t - tau)g(t) 가 더 이해하기 쉽다. g를 기준 신호로 보고, f에 tau만큼 딜레이를 줘가면서 유사도를 계산한 것이라 생각하면 편하기 때문.

wiener filter

https://dsp-nbsphinx.readthedocs.io/en/nbsphinx-experiment/randomsignalsLTIsystems/wienerfilter.html

• MSE로부터 wiener filter 최적해를 유도하는 과정은 건너뛰도록 하자…
• source의 distortion을 만드는 시스템 필터 G가 있을 경우, wiener filter H(w) = PSDsx(w) / PSDxx(w) 이다. 주파수 마다 계산된 다는 점에 주의.
• 이 PSDsx, PSDxx는 어떻게 아는가? 안다고 가정한다…

https://dsp-nbsphinx.readthedocs.io/en/nbsphinx-experiment/randomsignalsLTIsystems/wienerfilter.html#Transfer-Function-of-the-Wiener-Filter

• 여기 예제 코드에서 'shift for causal'이라는 주석으로 되어 있는 부분이 있는데, 이 간단한 개념을 이해하지 못해 시간 다 썼다.

H = H * np.exp(-1j*2*np.pi/len(H)*np.arange(len(H))*(len(H)//2)) # shift for causal

• '-1'는 반대 방향, 2pinp.arange(len(H))/len(H)는 frequency(w에 해당), (len(H)//2)는 shift해줄 길이(n또는 t)를 의미한다.

• 이는 코드상으로 구현한 wiener filter를 그대로 적용(convolution)할 경우 결과 신호가 반 파장 phase shift 되어 나타나는 것을 보정한다. non-causal filter라면 현재 위치 앞/뒤 신호를 한번에 보고 적용하면 되지만, causal system에 구현해야 할 경우 대칭의 중심이 필터 길이의 반에 위치하게 된다. causal filter로 얻은 결과 신호는 원 신호에 대해 lagging 되어 나타나기 때문에, 원 신호와 비교를 하기 위해서는 lagging 된 만큼 땡겨줄 필요가 있다.

• wiener deconvolution filter라 부르는 것은 wiener filter와 전혀 다를바가 없이 보이는데… 왜 이름을 다시 지었을까? 용도를 강조한 것 같다.

• wiener deconvolution filter는 wiener filter와 달리 PSDss, PSDnn, SNR, distortion filter G에 대한 수식으로 표현된다. 여기서 G(w)=1인 경우 H(w) = SNR(w) / (SNR(w) + 1) = PSDss / (PSDss + PSD_nn)이 된다.

• 가산 잡음이 없다면 H(w) = 1 / G(w) 가 된다.