配列の環状シフトなど

昨日「相互相関の計算 - 丸井綜研」で「circshiftのあたりが面倒な気がしたので」と書きましたが、そんなに面倒ではなかったので書き残しておきます。

(defun circshift (arry k)
  "Circular shift the array ARRY by amount K. Shift towards right if K
is positive number.
Example:
(circshift #(1 2 3 4 5)  1)
=> #(5 1 2 3 4)
(circshift #(1 2 3 4 5) -1)
;=> #(2 3 4 5 1)"
  (let ((res (make-array (length arry)))
        (amt (mod k (length arry))))
    (loop for n from amt to (1- (length arry))
          do (setf (aref res n) (aref arry (- n amt))))
    (loop for n from 0 to (1- amt)
          do (setf (aref res n) (aref arry (+ (- (length arry) amt) n))))
    res))

Matlabのcircshiftだと2次元以上の配列にも対応していますが、今回は簡易的に1次元配列だけに使えるようにしています。

これを使えば、昨日のMatlabコードは以下のように書くことができます。

(setq n (linspace 0 15 16))
(setq x (coerce (mapcar #'(lambda (x) (expt 0.84 x)) n) 'vector))
(setq y (circshift x 5))
(xcorr x y)

ついでに、似たような操作であるfftshiftとifftshiftも作りました。基本的には配列の左半分と右半分を入れ替えるだけなので、配列の長さが偶数であれば(fftshift (fftshift V))で元に戻ります。ただ、配列の長さが奇数の時には(fftshift (fftshift V))だとうまくいかないので、ifftshiftが必要になります。

(defun fftshift (arry)
  "Shift the FFT result to make the DC at the center of the array.
Example:
(fftshift #(1 2 3 4 5 6)
=> #(4 5 6 1 2 3)
(fftshift #(1 2 3 4 5 6 7)
=> #(5 6 7 1 2 3 4)"
  (if (evenp (length arry))
      (circshift arry (/ (length arry) 2))
      (circshift arry (/ (1- (length arry)) 2))))

(defun ifftshift (arry)
  "Inverse shift the FFT result to make the DC at the center of the array.
Example:
(ifftshift #(4 5 6 1 2 3)
=> #(1 2 3 4 5 6)
(ifftshift #(5 6 7 1 2 3 4)
=> #(1 2 3 4 5 6 7)"
  (if (evenp (length arry))
      (circshift arry (/ (length arry) 2))
      (circshift arry (/ (1+ (length arry)) 2))))