[SS] Circular Convolution

Circular Convolution이 필요한 이유.

cyclic convolution이라고도 불림.

• DFT의 경우, frequency domain representation을 샘플링하므로, time domain representation도 periodic signal이 됨.
• 때문에 DFT에서 time domain에서 input signal x[n]과 impulse response h[n]을 convolution하여 zero-state response y[n]를 구하는 경우, 기존의 (linear) convolution이 아닌 circular convolution으로 구해야 DFT의 sampling으로 인해 time domain의 x[n]과 h[n]이 주기를 가지고 반복되게 된 점을 반영할 수 있음.
• computer에서 신호를 다룬다는 건, 모든 신호가 discrete하게 존재한다는 것으로 결국 DFT와 같이 time domain이나 frequency domain이나 모두 periodic signal로 처리해야 한다.

DFT의 frequency domain에서의 sampling으로 인해 time domain의 signal(or function)들이 periodic이 되므로, linear convolution이 아닌 circular convolution이 이루어지게 됨.

Circular convolution 수행방법

길이가 4인 $x[n]$과 길이가 3인 $h[0]$를 circular convolution하면 다음과 같음.

1. x[n], h[n]의 크기를 같게 맞춤 (뒤에 0으로 padding)
   → 왼쪽 그림에서 h[3]=0으로 추가됨.
2. 왼쪽 그림(a)과 같이 원에 각 샘플 배치
3. 왼쪽 그림(b)와 같이 h[n]을 reflection (h[n]을 그대로 두고 x[n]을 시켜도 됨.)
4. 왼쪽 그림(c)와 같이 원의 중앙을 일치시키고 마주 보는 샘플들을 곱하고 그 결과들을 더해서 y[0]를 구함.
5. 이후 h[n]을 반시계방향으로 회전시키고, 4를 다시 수행하여 y[1]을 구함
6. y[n]의 샘플수가 1에서 맞춘 크기와 같아질 때까지 5를 수행.

참고사항

• 위 그림에서 h[n]이 오른쪽 테이블에서 이동하는 방식을 modular shift 또는 circular shift라고 함.
• (linear) convolution과 마찬가지로, circular convolution도 commutative임. (교환법칙 성립)
• 즉, 위 그림에서 h[n]과 x[n]의 처리를 바꾸어서 해도 상관없음.

$y[n]=x[n]\otimes h[n]$을 sliding table로 표현하면 다음과 같음.

Circular convolution과 DFT

DTFT의 convolution perporty와 마찬가지로, DFT에서도 time domain에서의 convolution은 frequency domain에서 곱에 해당함.

$$ y[n]=x[n]\otimes h[n] \underset{\text{DFT}}{\Longleftrightarrow} Y[k]=X[k]H[k] $$

때문에, x[n]과 h[n]에 대해 DFT를 수행하고, 결과 DFT들을 곱한 후 다시 IDFT를 하면 x[n]과 h[n]의 circular convolution y[n](=x[n]⊗h[n])을 구할 수 있음.

• 이 방법의 경우, 많은 라이브러리에서 최적화되어 제공하는 DFT, IDFT 수행 함수를 사용할 수 있음.

---

Python code examples

numpy 라이브러리를 이용하여 구현한다.

import numpy as np

크게 2가지로 구현을 함.

• ds_circular_conv1 : circular convolution의 알고리즘 그대로 구현하여 time domain에서 수행
• ds_circular_conv2 : DFT를 수행하여 frequenction domain에서 곱하기로 처리하는 방식을 이용한 것임.
• ds_circular_conv1은 circular_shift를 이용하는데, 알고리즘 그대로를 구현한 코드도 추가해둠.

하지만, 대부분의 경우 기존 라이브러리에서 해당 기능을 구현한 모듈이나 함수를 제공하며 이들을 이용하는게 성능이나 기능 면에서 항상 좋으므로 가급적 라이브러리에서 원하는 기능을 찾아보길 권함.

def circular_shift(signal_1d):
  """ circular shfit function for 1d signal
  example : [1 2 3 4 5] = [5 1 2 3 4]
  """

  """
  # 이해를 위한 구현.
  signal_size = len(signal_1d)
  ret = np.zeros(signal_size)

  last = signal_1d[-1]
  for i in range(1,signal_size):
    ret[i] = signal_1d[i-1]
    ret[0] = last
    return ret
  """
  # numpy 의 roll을 이용.  
  return np.roll(signal_1d,1)

  #  finding circular convolution
def ds_circular_conv1(x,h):
  """
  time domain에서 circular convolution 수행 
  (linear algebra의 matrix와 vector간 inner product를 이용해 구현함.)

  x : real input signal (1D real array)
  h : real impulse response (1D real array), kernel이라고도 불림.
  """
  if len(np.array(x).shape)>1 or len(np.array(x).shape)>1 :
    raise Exception('Invalid dimention error : Only 1-D data is available')
  input_size = len(x)
  res_size   = len(h)
  print(input_size, res_size)
  t_size = max(input_size, res_size)
  tmp = np.zeros( (t_size, t_size) )
  for i in range(input_size):
    tmp[0][i] = x[i]
  for i in range(1,t_size):
    tmp[i] = circular_shift(tmp[i-1])

  #print(tmp)

  transposed = tmp.T
  zero_padding = np.zeros((t_size-res_size))
  h_padded = np.concatenate( (np.array(h), zero_padding))
  ret = np.dot(transposed,h_padded.T)

  return ret

def ds_circular_conv2(x,h):
  """
  x : real input signal (1D real array)
  h : real impulse response (1D real array), kernel이라고도 불림.
  """
  if len(np.array(x).shape)>1 or len(np.array(x).shape)>1 :
    raise Exception('Invalid dimention error','Only 1-D data is available')

  input_size = len(x)
  res_size   = len(h)
  t_size = max(input_size, res_size)

  if(input_size > res_size):    
    h = np.pad(h, (0,t_size-res_size), 'constant') #, constant_values=(0,0))
  elif(input_size < res_size):
    x = np.pad(x, (0,t_size-input_size), 'constant') #, constant_values=(0,0))      

  return np.round(np.real(np.fft.ifft( np.fft.fft(x)*np.fft.fft(h) )),4)

---

Test : LPF

다음은 일종의 Low Pass Filter에 대한 결과를 확인한 것임.

x[n] 의 길이가 8이고 h[n] 의 길이가 4임.
x[n] 의 주기 N=2 임.
Low pass filter로 DC성분의 크기 5로 출력신호 (DC에 대해 gain이 1임)가 나오게 됨.

x = [10,0,10,0,10,0,10,0]
h = [1/4,1/4,1/4,1/4]

circular_convolution_result = ds_circular_conv1(x,h)
print(circular_convolution_result)

결과는 다음과 같이 DC성분으로 나옴.

[5. 5. 5. 5. 5. 5. 5. 5.]

---

Test : HPF

다음은 일종의 High Pass Filter에 대한 결과를 확인한 것임.

• x[n] 의 길이가 8이고 h[n] 의 길이가 4임.
• x[n] 의 주기 N=2 임.
  High pass filter로 DC성분인 5로 일정한 부분이 빠지고 amplitude 5로 변화하는 고주파 성분(최고치와 최저치의 차이는 5*2=10)이 출력신호로 나옴

결과는 다음과 같이 DC가 빠져서 고주파 성분만이 나옴.
[ 5. -5. 5. -5. 5. -5. 5. -5.]
source code는 다음과 같음.

x = [10,0,10,0,10,0,10,0]
h = [1/4,-1/4,1/4,-1/4]
circular_convolution_result = ds_circular_conv1(x,h) 
print(circular_convolution_result)

---

Circular convolution으로 (linear) convolution과 같은 결과 얻기.

일반적으로 linear convolution과 cyclic convolution은 결과가 같지 않으므로, circular convolution으로 linear convolution과 같은 결과를 얻으려면 operand에 해당하는 $x[n]$과 $h[n]$을 적절히 zero padding해야 한다.

구체적인 알고리즘은 다음과 같음.

• 우선 convolution의 operand인 $x[n]$과 $h[n]$을 l=len(x)+len(h)-1의 크기로 zero padding을 수행함.
• 이후 circular convolution을 수행하면 (linear) convolution의 full mode 결과가 나옴.
• valid나 same은 각 모드에 따라 full mode 결과에서 앞뒤를 잘라서 사용하면 됨.

code로 확인하면 다음과 같음.

x = [10,0,10,0,10,0,10,0]
h = [1/4,-1/4,1/4,-1/4]

print(np.convolve(x,h, mode='full')) # linear convolution

l = len(x)+len(h)-1

x_p = np.pad(x,(0,l-len(x)),mode='constant')
h_p = np.pad(h,(0,l-len(h)),mode='constant')
print(ds_circular_conv2(x_p,h_p)) # circular convolution

---

관련 ipynb

https://gist.github.com/dsaint31x/47caa2064184bb925587bb05888bc6e1

Circular_Convolution_N_Linear_Convolution.ipynb