오버플로 없는 short int x 분수

이슈

저의 OAGLE에는 오디오 파일을 복호화해서 재생하는 부분이 있습니다. 단순히 재생하는 것 외에 대략적으로는 다음과 같은 기능적 요구 사항이 있습니다.

• 여러 개의 오디오 파일을 동시에 재생한다
• 각 재생 스트림의 볼륨을 따로 조정한다
• 마스터 볼륨

단순히 하나의 파일을 재생한다면 memcpy()의 최적화가 아주 훌륭하므로 복호화한 걸 그대로 링버퍼에 갖다가 붙이면 그만입니다. 전자는 PCM끼리 더하면 되고 후자는 PCM에 볼륨을 곱해야 합니다. 이 과정은 memcpy()에 비해 오래 걸립니다. 초당 44100 샘플이 필요하다면 SSE2 float 기준으로 초당 11025 x (n-1)회의 실수배와 벡터 합 연산이 필요한 셈입니다. 그럼에도 저는 프로그램 상의 샘플 기준으로 float를 쓰고 있었는데요, 오버플로 관리가 필요 없기 때문입니다. float 샘플은 [-1,1] 범위가 기준이 되니 다 더하고 마지막에 1로 컷하면 된다는 겁니다.

short 타입의 샘플은 short 범위 [-32768, 32767]를 모두 사용합니다. SSE2에서는 _mm_adds_epi16 함수가 있고 이는 일반 합과 성능 차이가 거의 없지만 범위를 제한해 줍니다(saturated arithmetic).

일단 인텔의 문서에 따르면, 필요한 요소의 처리율(Skylake 기준)은 다음과 같습니다. 처리율의 단위는 CPI(파이프라인 상에서 명령당 사이클)로, 이 값은 낮을수록 좋습니다(오타 아님). 사실 표를 만들 필요도 없이 사용할 intrinsic 함수의 이론상 처리율은 같습니다.

intrinsic 함수	CPI
_mm_add_ps	0.5
_mm_mul_ps	0.5
_mm_mulhi_epi16 _mm_mullo_epi16	0.5
_mm_adds_epi16	0.5
_mm_srai_epi16	0.5
_mm_loadu_ps _mm_loadu_si128	0.5

 

float로 처리할 때, 볼륨을 곱해서 더해지는 스트림마다 처리할 PCM 샘플이 8a개라고 할 때 성능은 이렇습니다.

- 곱 2a회

- 합 2a회

- 로드 4a+1회 (상수 로드 속도를 같다고 가정)

총 8a+1

 

볼륨이 0 혹은 100%라면 곱할 이유가 없습니다. 이는 분기를 나눠서라도 피해야 합니다.

그냥 더해지는 스트림에 대한 성능은 이렇습니다.

- 합 2a회

- 로드 4a회

총 6a

 

그냥 더해지는 스트림을 short로 처리하면 성능은 다음과 같습니다.

- 합 a회

- 로드 2a회

총 3a

 

볼륨을 1로 둔다면 2배 이익이니, 볼륨을 낮게 조정했을 때의 필요 사이클 단위가 8a를 넘지만 않으면 충분히 바꿀 가치가 있겠습니다. 하지만 문제는 제목에서와 같습니다. 여러 개 재생은 그렇다 치고 볼륨 조절은 어떻게 할까요? short를 32비트 혹은 64비트 자료형으로, 그리고 그 반대로 바꾸는 것은 오버헤드가 없는 수준이긴 합니다만 SIMD가 들어가면 말이 달라집니다. 확장하고 나서 더한다면 float를 그대로 쓰는 것에 비해 이득이 크지 않죠. 외려 32비트로 올라간 값을 다시 16비트 오프셋에 맞춰 담는 데 시간을 쓸 것 같습니다. 그러니 어떤 방법을 써서든 이 볼륨 곱하기를 구현했으면 float와 성능을 비교했을 때 좀 괜찮았으면 좋겠습니다.

문제: 정수 x 분수를 최대한 적은 오차로 계산하기

당연히 부동소수점 타입으로 변환하거나 union으로 겹치는 것도 32비트를 써야 하므로 논외입니다. 물론 그렇대도 답은 간단합니다. 중점은 이론상 계산량 비교가 되겠습니다.

일반적인 분수를 다 쓴다면 무슨 일이 있어도 float보다 느릴 수밖에 없으니 분모는 2ⁿ으로 고정해 줍시다. 볼륨은 보통 백분율을 사용하므로 128 이상의 분모는 대체로 충분한 정밀도를 제공합니다.

1. 1/2ⁿ 곱하기

이건 >> (산술)시프트 한 번으로 충분합니다. _mm_srai_epi16이 그에 대응하며, 오차는 당연히 1 미만입니다. 이것은 문제의 직접적 해결은 아니므로 성능에 대한 이야기는 의미가 없어서 하지 않겠습니다.

2. m/2ⁿ 곱하기 (m < 2ⁿ)

만약 산술 시프트 한 번 후 m을 곱하는 것으로 끝낸다면, 오차는 m 미만이 되겠죠. 마찬가지로 처리할 PCM 샘플이 8a개라고 할 때 성능은 이렇습니다.

- 시프트 a회

- 곱 a회

- 포화합 a회

- 로드 2a+1회

총 5a+1

나쁘지 않습니다. 하지만 오차 범위는 꽤 나쁩니다. 분모가 128이라면 오차는 최대 127. (float로 치면 0.0038 정도)

 

오차가 없는 것 중 가장 쉬운 방법은 분모를 65536으로 고정하고, mulhi를 쓰는 겁니다. 이 경우 이미 우측으로 16 시프트한 것이라는 차이가 있습니다. 코드로 치면 아래와 같습니다.

__m128i sh = _mm_set1_epi16(32000);	// 32000/65536 을 곱하기
xx = _mm_mulhi_epi16(yy, sh);

이 경우,

- 곱 a회

- 포화합 a회

- 로드 2a+1회

총 4a+1이라는, 오차가 1 미만인 주제에 위보다 강력한 성능이 나옵니다. 하지만 위에는 없는 문제가 여기에는 있습니다. 곱셈의 결과가 제대로 나올 수 있는 것은 +32767까지라는 겁니다. 이것을 해결하는 방법은 2가지가 생각납니다.

 

1. 분모를 65536 말고 32768로 하며, 이에 따라 결과에 좌측 시프트를 한 번 추가로 한다(2를 곱한다)

이 경우 시프트를 또 해야 하니 5a+1의 성능을 가지게 되며 오차는 최대 2가 됩니다. 2 정도면 거의 무시할 수 있는 수준이 되긴 합니다.

 

2. 32767이 넘으면 빼기로 전환한다

볼륨값이 정해졌을 때 1/2이 넘는다면 원본에서 빼기로 전환하는 방법입니다. 1/2을 넘는 경우 역시 성능은 5a+1입니다. 예를 들어 32767에 short인 65000(실제 값 -536)을 곱한 상위 16비트의 결과는 -268로, 이것을 원본에 더하면 32499입니다. 실제 32767 x 65000 / 65536 = 32499.0082 정도입니다.

스트림당 분기가 한 번이므로, 분기는 성능에 영향을 거의 미치지 않습니다. 아마 이 말이 이해가 잘 안 될 수도 있는데, 링 버퍼를 구현해 본 사람을 위한 설명입니다. 대단한 건 아니고 볼륨이 정해진 순간 1/2과 비교하여 그 스트림에 대한 플래그가 값 65000과 함께 저장된다는 말입니다.

 

* 부호형과 비부호형 정수의 계산 상 취급에는 원래 차이가 없습니다. 하지만 크기가 확장되는 명령에서는 예외인데요, 예를 들어 -32768은 비트가 1000 0000 0000 0000인데 비부호형 계산에서는 0000 0000 0000 0000 1000 0000 0000 0000으로 확장되지만 부호형 계산에서는 1111 1111 1111 1111 1000 0000 0000 0000으로 확장됩니다. 그래서 곱해지는 수만 비부호형으로 취급하고 그런 일이 불가능합니다.

 

속도로 치면 아마 2번이 조금 더 빠를 것 같습니다. 하지만 코드가 살짝 지저분해질 것 같긴 해요.