diff --git a/av1/av1.cmake b/av1/av1.cmake
index b84b83b..7976746 100644
--- a/av1/av1.cmake
+++ b/av1/av1.cmake
@@ -274,6 +274,7 @@
             "${AOM_ROOT}/av1/encoder/x86/corner_match_sse4.c"
             "${AOM_ROOT}/av1/encoder/x86/encodetxb_sse4.c"
             "${AOM_ROOT}/av1/encoder/x86/highbd_fwd_txfm_sse4.c"
+            "${AOM_ROOT}/av1/encoder/x86/rdopt_sse4.c"
             "${AOM_ROOT}/av1/encoder/x86/pickrst_sse4.c")
 
 list(APPEND AOM_AV1_ENCODER_INTRIN_AVX2
@@ -284,6 +285,7 @@
             "${AOM_ROOT}/av1/encoder/x86/av1_fwd_txfm2d_avx2.c"
             "${AOM_ROOT}/av1/encoder/x86/wedge_utils_avx2.c"
             "${AOM_ROOT}/av1/encoder/x86/encodetxb_avx2.c"
+            "${AOM_ROOT}/av1/encoder/x86/rdopt_avx2.c"
             "${AOM_ROOT}/av1/encoder/x86/pickrst_avx2.c")
 
 list(APPEND AOM_AV1_ENCODER_INTRIN_NEON
diff --git a/av1/common/av1_rtcd_defs.pl b/av1/common/av1_rtcd_defs.pl
index aabb001..f2a043c 100755
--- a/av1/common/av1_rtcd_defs.pl
+++ b/av1/common/av1_rtcd_defs.pl
@@ -305,6 +305,9 @@
 
   add_proto qw/int64_t av1_highbd_pixel_proj_error/, " const uint8_t *src8, int width, int height, int src_stride, const uint8_t *dat8, int dat_stride, int32_t *flt0, int flt0_stride, int32_t *flt1, int flt1_stride, int xq[2], const sgr_params_type *params";
   specialize qw/av1_highbd_pixel_proj_error sse4_1 avx2/;
+
+  add_proto qw/void av1_get_horver_correlation_full/, " const int16_t *diff, int stride, int w, int h, float *hcorr, float *vcorr";
+  specialize qw/av1_get_horver_correlation_full sse4_1 avx2/;
 }
 # end encoder functions
 
diff --git a/av1/encoder/rdopt.c b/av1/encoder/rdopt.c
index 970260e..e27a713 100644
--- a/av1/encoder/rdopt.c
+++ b/av1/encoder/rdopt.c
@@ -1669,69 +1669,100 @@
 
 // Similar to get_horver_correlation, but also takes into account first
 // row/column, when computing horizontal/vertical correlation.
-static void get_horver_correlation_full(const int16_t *diff, int stride, int w,
-                                        int h, float *hcorr, float *vcorr) {
-  const float num_hor = (float)(h * (w - 1));
-  const float num_ver = (float)((h - 1) * w);
-  int i, j;
-
+void av1_get_horver_correlation_full_c(const int16_t *diff, int stride,
+                                       int width, int height, float *hcorr,
+                                       float *vcorr) {
   // The following notation is used:
   // x - current pixel
   // y - left neighbor pixel
   // z - top neighbor pixel
-  int64_t xy_sum = 0, xz_sum = 0;
-  int64_t xhor_sum = 0, xver_sum = 0, y_sum = 0, z_sum = 0;
-  int64_t x2hor_sum = 0, x2ver_sum = 0, y2_sum = 0, z2_sum = 0;
+  int64_t x_sum = 0, x2_sum = 0, xy_sum = 0, xz_sum = 0;
+  int64_t x_firstrow = 0, x_finalrow = 0, x_firstcol = 0, x_finalcol = 0;
+  int64_t x2_firstrow = 0, x2_finalrow = 0, x2_firstcol = 0, x2_finalcol = 0;
 
-  int16_t x, y, z;
-  for (j = 1; j < w; ++j) {
-    x = diff[j];
-    y = diff[j - 1];
+  // First, process horizontal correlation on just the first row
+  x_sum += diff[0];
+  x2_sum += diff[0] * diff[0];
+  x_firstrow += diff[0];
+  x2_firstrow += diff[0] * diff[0];
+  for (int j = 1; j < width; ++j) {
+    const int16_t x = diff[j];
+    const int16_t y = diff[j - 1];
+    x_sum += x;
+    x_firstrow += x;
+    x2_sum += x * x;
+    x2_firstrow += x * x;
     xy_sum += x * y;
-    xhor_sum += x;
-    y_sum += y;
-    x2hor_sum += x * x;
-    y2_sum += y * y;
   }
-  for (i = 1; i < h; ++i) {
-    x = diff[i * stride];
-    z = diff[(i - 1) * stride];
+
+  // Process vertical correlation in the first column
+  x_firstcol += diff[0];
+  x2_firstcol += diff[0] * diff[0];
+  for (int i = 1; i < height; ++i) {
+    const int16_t x = diff[i * stride];
+    const int16_t z = diff[(i - 1) * stride];
+    x_sum += x;
+    x_firstcol += x;
+    x2_sum += x * x;
+    x2_firstcol += x * x;
     xz_sum += x * z;
-    xver_sum += x;
-    z_sum += z;
-    x2ver_sum += x * x;
-    z2_sum += z * z;
-    for (j = 1; j < w; ++j) {
-      x = diff[i * stride + j];
-      y = diff[i * stride + j - 1];
-      z = diff[(i - 1) * stride + j];
+  }
+
+  // Now process horiz and vert correlation through the rest unit
+  for (int i = 1; i < height; ++i) {
+    for (int j = 1; j < width; ++j) {
+      const int16_t x = diff[i * stride + j];
+      const int16_t y = diff[i * stride + j - 1];
+      const int16_t z = diff[(i - 1) * stride + j];
+      x_sum += x;
+      x2_sum += x * x;
       xy_sum += x * y;
       xz_sum += x * z;
-      xhor_sum += x;
-      xver_sum += x;
-      y_sum += y;
-      z_sum += z;
-      x2hor_sum += x * x;
-      x2ver_sum += x * x;
-      y2_sum += y * y;
-      z2_sum += z * z;
     }
   }
+
+  for (int j = 0; j < width; ++j) {
+    x_finalrow += diff[(height - 1) * stride + j];
+    x2_finalrow +=
+        diff[(height - 1) * stride + j] * diff[(height - 1) * stride + j];
+  }
+  for (int i = 0; i < height; ++i) {
+    x_finalcol += diff[i * stride + width - 1];
+    x2_finalcol += diff[i * stride + width - 1] * diff[i * stride + width - 1];
+  }
+
+  int64_t xhor_sum = x_sum - x_finalcol;
+  int64_t xver_sum = x_sum - x_finalrow;
+  int64_t y_sum = x_sum - x_firstcol;
+  int64_t z_sum = x_sum - x_firstrow;
+  int64_t x2hor_sum = x2_sum - x2_finalcol;
+  int64_t x2ver_sum = x2_sum - x2_finalrow;
+  int64_t y2_sum = x2_sum - x2_firstcol;
+  int64_t z2_sum = x2_sum - x2_firstrow;
+
+  const float num_hor = (float)(height * (width - 1));
+  const float num_ver = (float)((height - 1) * width);
+
   const float xhor_var_n = x2hor_sum - (xhor_sum * xhor_sum) / num_hor;
-  const float y_var_n = y2_sum - (y_sum * y_sum) / num_hor;
-  const float xy_var_n = xy_sum - (xhor_sum * y_sum) / num_hor;
   const float xver_var_n = x2ver_sum - (xver_sum * xver_sum) / num_ver;
+
+  const float y_var_n = y2_sum - (y_sum * y_sum) / num_hor;
   const float z_var_n = z2_sum - (z_sum * z_sum) / num_ver;
+
+  const float xy_var_n = xy_sum - (xhor_sum * y_sum) / num_hor;
   const float xz_var_n = xz_sum - (xver_sum * z_sum) / num_ver;
 
-  *hcorr = *vcorr = 1;
   if (xhor_var_n > 0 && y_var_n > 0) {
     *hcorr = xy_var_n / sqrtf(xhor_var_n * y_var_n);
     *hcorr = *hcorr < 0 ? 0 : *hcorr;
+  } else {
+    *hcorr = 1.0;
   }
   if (xver_var_n > 0 && z_var_n > 0) {
     *vcorr = xz_var_n / sqrtf(xver_var_n * z_var_n);
     *vcorr = *vcorr < 0 ? 0 : *vcorr;
+  } else {
+    *vcorr = 1.0;
   }
 }
 
@@ -1844,9 +1875,9 @@
   const int16_t *diff = p->src_diff + 4 * blk_row * diff_stride + 4 * blk_col;
   get_energy_distribution_finer(diff, diff_stride, bw, bh, hfeatures,
                                 vfeatures);
-  get_horver_correlation_full(diff, diff_stride, bw, bh,
-                              &hfeatures[hfeatures_num - 1],
-                              &vfeatures[vfeatures_num - 1]);
+  av1_get_horver_correlation_full(diff, diff_stride, bw, bh,
+                                  &hfeatures[hfeatures_num - 1],
+                                  &vfeatures[vfeatures_num - 1]);
   av1_nn_predict(hfeatures, nn_config_hor, hscores);
   av1_nn_predict(vfeatures, nn_config_ver, vscores);
 
diff --git a/av1/encoder/x86/rdopt_avx2.c b/av1/encoder/x86/rdopt_avx2.c
new file mode 100644
index 0000000..a94c076
--- /dev/null
+++ b/av1/encoder/x86/rdopt_avx2.c
@@ -0,0 +1,253 @@
+/*
+ * Copyright (c) 2018, Alliance for Open Media. All rights reserved
+ *
+ * This source code is subject to the terms of the BSD 2 Clause License and
+ * the Alliance for Open Media Patent License 1.0. If the BSD 2 Clause License
+ * was not distributed with this source code in the LICENSE file, you can
+ * obtain it at www.aomedia.org/license/software. If the Alliance for Open
+ * Media Patent License 1.0 was not distributed with this source code in the
+ * PATENTS file, you can obtain it at www.aomedia.org/license/patent.
+ */
+
+#include <assert.h>
+#include <immintrin.h>
+#include "aom_dsp/x86/synonyms_avx2.h"
+
+#include "config/av1_rtcd.h"
+#include "av1/encoder/rdopt.h"
+
+// Process horizontal and vertical correlations in a 4x4 block of pixels.
+// We actually use the 4x4 pixels to calculate correlations corresponding to
+// the top-left 3x3 pixels, so this function must be called with 1x1 overlap,
+// moving the window along/down by 3 pixels at a time.
+INLINE static void horver_correlation_4x4(const int16_t *diff, int stride,
+                                          __m256i *xy_sum_32,
+                                          __m256i *xz_sum_32, __m256i *x_sum_32,
+                                          __m256i *x2_sum_32) {
+  // Pixels in this 4x4   [ a b c d ]
+  // are referred to as:  [ e f g h ]
+  //                      [ i j k l ]
+  //                      [ m n o p ]
+
+  const __m256i pixels = _mm256_set_epi64x(
+      *(uint64_t *)&diff[0 * stride], *(uint64_t *)&diff[1 * stride],
+      *(uint64_t *)&diff[2 * stride], *(uint64_t *)&diff[3 * stride]);
+  // pixels = [d c b a h g f e] [l k j i p o n m] as i16
+
+  const __m256i slli = _mm256_slli_epi64(pixels, 16);
+  // slli = [c b a 0 g f e 0] [k j i 0 o n m 0] as i16
+
+  const __m256i madd_xy = _mm256_madd_epi16(pixels, slli);
+  // madd_xy = [bc+cd ab fg+gh ef] [jk+kl ij no+op mn] as i32
+  *xy_sum_32 = _mm256_add_epi32(*xy_sum_32, madd_xy);
+
+  // Permute control [3 2] [1 0] => [2 1] [0 0], 0b10010000 = 0x90
+  const __m256i perm = _mm256_permute4x64_epi64(slli, 0x90);
+  // perm = [g f e 0 k j i 0] [o n m 0 o n m 0] as i16
+
+  const __m256i madd_xz = _mm256_madd_epi16(slli, perm);
+  // madd_xz = [cg+bf ae gk+fj ei] [ko+jn im oo+nn mm] as i32
+  *xz_sum_32 = _mm256_add_epi32(*xz_sum_32, madd_xz);
+
+  // Sum every element in slli (and then also their squares)
+  const __m256i madd1_slli = _mm256_madd_epi16(slli, _mm256_set1_epi16(1));
+  // madd1_slli = [c+b a g+f e] [k+j i o+n m] as i32
+  *x_sum_32 = _mm256_add_epi32(*x_sum_32, madd1_slli);
+
+  const __m256i madd_slli = _mm256_madd_epi16(slli, slli);
+  // madd_slli = [cc+bb aa gg+ff ee] [kk+jj ii oo+nn mm] as i32
+  *x2_sum_32 = _mm256_add_epi32(*x2_sum_32, madd_slli);
+}
+
+void av1_get_horver_correlation_full_avx2(const int16_t *diff, int stride,
+                                          int width, int height, float *hcorr,
+                                          float *vcorr) {
+  // The following notation is used:
+  // x - current pixel
+  // y - right neighbour pixel
+  // z - below neighbour pixel
+  // w - down-right neighbour pixel
+  int64_t xy_sum = 0, xz_sum = 0;
+  int64_t x_sum = 0, x2_sum = 0;
+
+  // Process horizontal and vertical correlations through the body in 4x4
+  // blocks.  This excludes the final row and column and possibly one extra
+  // column depending how 3 divides into width and height
+  int32_t xy_xz_tmp[8] = { 0 }, x_x2_tmp[8] = { 0 };
+  __m256i xy_sum_32 = _mm256_setzero_si256();
+  __m256i xz_sum_32 = _mm256_setzero_si256();
+  __m256i x_sum_32 = _mm256_setzero_si256();
+  __m256i x2_sum_32 = _mm256_setzero_si256();
+  for (int i = 0; i <= height - 4; i += 3) {
+    for (int j = 0; j <= width - 4; j += 3) {
+      horver_correlation_4x4(&diff[i * stride + j], stride, &xy_sum_32,
+                             &xz_sum_32, &x_sum_32, &x2_sum_32);
+    }
+    const __m256i hadd_xy_xz = _mm256_hadd_epi32(xy_sum_32, xz_sum_32);
+    // hadd_xy_xz = [ae+bf+cg ei+fj+gk ab+bc+cd ef+fg+gh]
+    //              [im+jn+ko mm+nn+oo ij+jk+kl mn+no+op] as i32
+    yy_storeu_256(xy_xz_tmp, hadd_xy_xz);
+    xy_sum += (int64_t)xy_xz_tmp[5] + xy_xz_tmp[4] + xy_xz_tmp[1];
+    xz_sum += (int64_t)xy_xz_tmp[7] + xy_xz_tmp[6] + xy_xz_tmp[3];
+
+    const __m256i hadd_x_x2 = _mm256_hadd_epi32(x_sum_32, x2_sum_32);
+    // hadd_x_x2 = [aa+bb+cc ee+ff+gg a+b+c e+f+g]
+    //             [ii+jj+kk mm+nn+oo i+j+k m+n+o] as i32
+    yy_storeu_256(x_x2_tmp, hadd_x_x2);
+    x_sum += (int64_t)x_x2_tmp[5] + x_x2_tmp[4] + x_x2_tmp[1];
+    x2_sum += (int64_t)x_x2_tmp[7] + x_x2_tmp[6] + x_x2_tmp[3];
+
+    xy_sum_32 = _mm256_setzero_si256();
+    xz_sum_32 = _mm256_setzero_si256();
+    x_sum_32 = _mm256_setzero_si256();
+    x2_sum_32 = _mm256_setzero_si256();
+  }
+
+  // x_sum now covers every pixel except the final 1-2 rows and 1-2 cols
+  int64_t x_finalrow = 0, x_finalcol = 0, x2_finalrow = 0, x2_finalcol = 0;
+
+  // Do we have 2 rows remaining or just the one?  Note that width and height
+  // are powers of 2, so each modulo 3 must be 1 or 2.
+  if (height % 3 == 1) {  // Just horiz corrs on the final row
+    const int16_t x0 = diff[(height - 1) * stride];
+    x_sum += x0;
+    x_finalrow += x0;
+    x2_sum += x0 * x0;
+    x2_finalrow += x0 * x0;
+    for (int j = 0; j < width - 1; ++j) {
+      const int16_t x = diff[(height - 1) * stride + j];
+      const int16_t y = diff[(height - 1) * stride + j + 1];
+      xy_sum += x * y;
+      x_sum += y;
+      x2_sum += y * y;
+      x_finalrow += y;
+      x2_finalrow += y * y;
+    }
+  } else {  // Two rows remaining to do
+    const int16_t x0 = diff[(height - 2) * stride];
+    const int16_t z0 = diff[(height - 1) * stride];
+    x_sum += x0 + z0;
+    x2_sum += x0 * x0 + z0 * z0;
+    x_finalrow += z0;
+    x2_finalrow += z0 * z0;
+    for (int j = 0; j < width - 1; ++j) {
+      const int16_t x = diff[(height - 2) * stride + j];
+      const int16_t y = diff[(height - 2) * stride + j + 1];
+      const int16_t z = diff[(height - 1) * stride + j];
+      const int16_t w = diff[(height - 1) * stride + j + 1];
+
+      // Horizontal and vertical correlations for the penultimate row:
+      xy_sum += x * y;
+      xz_sum += x * z;
+
+      // Now just horizontal correlations for the final row:
+      xy_sum += z * w;
+
+      x_sum += y + w;
+      x2_sum += y * y + w * w;
+      x_finalrow += w;
+      x2_finalrow += w * w;
+    }
+  }
+
+  // Do we have 2 columns remaining or just the one?
+  if (width % 3 == 1) {  // Just vert corrs on the final col
+    const int16_t x0 = diff[width - 1];
+    x_sum += x0;
+    x_finalcol += x0;
+    x2_sum += x0 * x0;
+    x2_finalcol += x0 * x0;
+    for (int i = 0; i < height - 1; ++i) {
+      const int16_t x = diff[i * stride + width - 1];
+      const int16_t z = diff[(i + 1) * stride + width - 1];
+      xz_sum += x * z;
+      x_finalcol += z;
+      x2_finalcol += z * z;
+      // So the bottom-right elements don't get counted twice:
+      if (i < height - (height % 3 == 1 ? 2 : 3)) {
+        x_sum += z;
+        x2_sum += z * z;
+      }
+    }
+  } else {  // Two cols remaining
+    const int16_t x0 = diff[width - 2];
+    const int16_t y0 = diff[width - 1];
+    x_sum += x0 + y0;
+    x2_sum += x0 * x0 + y0 * y0;
+    x_finalcol += y0;
+    x2_finalcol += y0 * y0;
+    for (int i = 0; i < height - 1; ++i) {
+      const int16_t x = diff[i * stride + width - 2];
+      const int16_t y = diff[i * stride + width - 1];
+      const int16_t z = diff[(i + 1) * stride + width - 2];
+      const int16_t w = diff[(i + 1) * stride + width - 1];
+
+      // Horizontal and vertical correlations for the penultimate col:
+      // Skip these on the last iteration of this loop if we also had two
+      // rows remaining, otherwise the final horizontal and vertical correlation
+      // get erroneously processed twice
+      if (i < height - 2 || height % 3 == 1) {
+        xy_sum += x * y;
+        xz_sum += x * z;
+      }
+
+      x_finalcol += w;
+      x2_finalcol += w * w;
+      // So the bottom-right elements don't get counted twice:
+      if (i < height - (height % 3 == 1 ? 2 : 3)) {
+        x_sum += z + w;
+        x2_sum += z * z + w * w;
+      }
+
+      // Now just vertical correlations for the final column:
+      xz_sum += y * w;
+    }
+  }
+
+  // Calculate the simple sums and squared-sums
+  int64_t x_firstrow = 0, x_firstcol = 0;
+  int64_t x2_firstrow = 0, x2_firstcol = 0;
+
+  for (int j = 0; j < width; ++j) {
+    x_firstrow += diff[j];
+    x2_firstrow += diff[j] * diff[j];
+  }
+  for (int i = 0; i < height; ++i) {
+    x_firstcol += diff[i * stride];
+    x2_firstcol += diff[i * stride] * diff[i * stride];
+  }
+
+  int64_t xhor_sum = x_sum - x_finalcol;
+  int64_t xver_sum = x_sum - x_finalrow;
+  int64_t y_sum = x_sum - x_firstcol;
+  int64_t z_sum = x_sum - x_firstrow;
+  int64_t x2hor_sum = x2_sum - x2_finalcol;
+  int64_t x2ver_sum = x2_sum - x2_finalrow;
+  int64_t y2_sum = x2_sum - x2_firstcol;
+  int64_t z2_sum = x2_sum - x2_firstrow;
+
+  const float num_hor = (float)(height * (width - 1));
+  const float num_ver = (float)((height - 1) * width);
+
+  const float xhor_var_n = x2hor_sum - (xhor_sum * xhor_sum) / num_hor;
+  const float xver_var_n = x2ver_sum - (xver_sum * xver_sum) / num_ver;
+
+  const float y_var_n = y2_sum - (y_sum * y_sum) / num_hor;
+  const float z_var_n = z2_sum - (z_sum * z_sum) / num_ver;
+
+  const float xy_var_n = xy_sum - (xhor_sum * y_sum) / num_hor;
+  const float xz_var_n = xz_sum - (xver_sum * z_sum) / num_ver;
+
+  if (xhor_var_n > 0 && y_var_n > 0) {
+    *hcorr = xy_var_n / sqrtf(xhor_var_n * y_var_n);
+    *hcorr = *hcorr < 0 ? 0 : *hcorr;
+  } else {
+    *hcorr = 1.0;
+  }
+  if (xver_var_n > 0 && z_var_n > 0) {
+    *vcorr = xz_var_n / sqrtf(xver_var_n * z_var_n);
+    *vcorr = *vcorr < 0 ? 0 : *vcorr;
+  } else {
+    *vcorr = 1.0;
+  }
+}
diff --git a/av1/encoder/x86/rdopt_sse4.c b/av1/encoder/x86/rdopt_sse4.c
new file mode 100644
index 0000000..f5ffae7
--- /dev/null
+++ b/av1/encoder/x86/rdopt_sse4.c
@@ -0,0 +1,272 @@
+/*
+ * Copyright (c) 2018, Alliance for Open Media. All rights reserved
+ *
+ * This source code is subject to the terms of the BSD 2 Clause License and
+ * the Alliance for Open Media Patent License 1.0. If the BSD 2 Clause License
+ * was not distributed with this source code in the LICENSE file, you can
+ * obtain it at www.aomedia.org/license/software. If the Alliance for Open
+ * Media Patent License 1.0 was not distributed with this source code in the
+ * PATENTS file, you can obtain it at www.aomedia.org/license/patent.
+ */
+
+#include <assert.h>
+#include <emmintrin.h>
+#include "aom_dsp/x86/synonyms.h"
+
+#include "config/av1_rtcd.h"
+#include "av1/encoder/rdopt.h"
+
+// Process horizontal and vertical correlations in a 4x4 block of pixels.
+// We actually use the 4x4 pixels to calculate correlations corresponding to
+// the top-left 3x3 pixels, so this function must be called with 1x1 overlap,
+// moving the window along/down by 3 pixels at a time.
+INLINE static void horver_correlation_4x4(const int16_t *diff, int stride,
+                                          __m128i *xy_sum_32,
+                                          __m128i *xz_sum_32, __m128i *x_sum_32,
+                                          __m128i *x2_sum_32) {
+  // Pixels in this 4x4   [ a b c d ]
+  // are referred to as:  [ e f g h ]
+  //                      [ i j k l ]
+  //                      [ m n o p ]
+
+  const __m128i pixelsa = _mm_set_epi64x(*(uint64_t *)&diff[0 * stride],
+                                         *(uint64_t *)&diff[2 * stride]);
+  const __m128i pixelsb = _mm_set_epi64x(*(uint64_t *)&diff[1 * stride],
+                                         *(uint64_t *)&diff[3 * stride]);
+  // pixelsa = [d c b a l k j i] as i16
+  // pixelsb = [h g f e p o n m] as i16
+
+  const __m128i slli_a = _mm_slli_epi64(pixelsa, 16);
+  const __m128i slli_b = _mm_slli_epi64(pixelsb, 16);
+  // slli_a = [c b a 0 k j i 0] as i16
+  // slli_b = [g f e 0 o n m 0] as i16
+
+  const __m128i xy_madd_a = _mm_madd_epi16(pixelsa, slli_a);
+  const __m128i xy_madd_b = _mm_madd_epi16(pixelsb, slli_b);
+  // xy_madd_a = [bc+cd ab jk+kl ij] as i32
+  // xy_madd_b = [fg+gh ef no+op mn] as i32
+
+  const __m128i xy32 = _mm_hadd_epi32(xy_madd_b, xy_madd_a);
+  // xy32 = [ab+bc+cd ij+jk+kl ef+fg+gh mn+no+op] as i32
+  *xy_sum_32 = _mm_add_epi32(*xy_sum_32, xy32);
+
+  const __m128i xz_madd_a = _mm_madd_epi16(slli_a, slli_b);
+  // xz_madd_a = [bf+cg ae jn+ko im] i32
+
+  const __m128i swap_b = _mm_srli_si128(slli_b, 8);
+  // swap_b = [0 0 0 0 g f e 0] as i16
+  const __m128i xz_madd_b = _mm_madd_epi16(slli_a, swap_b);
+  // xz_madd_b = [0 0 gk+fj ei] i32
+
+  const __m128i xz32 = _mm_hadd_epi32(xz_madd_b, xz_madd_a);
+  // xz32 = [ae+bf+cg im+jn+ko 0 ei+fj+gk] i32
+  *xz_sum_32 = _mm_add_epi32(*xz_sum_32, xz32);
+
+  // Now calculate the straight sums, x_sum += a+b+c+e+f+g+i+j+k
+  // (sum up every element in slli_a and swap_b)
+  const __m128i sum_slli_a = _mm_hadd_epi16(slli_a, slli_a);
+  const __m128i sum_slli_a32 = _mm_cvtepi16_epi32(sum_slli_a);
+  // sum_slli_a32 = [c+b a k+j i] as i32
+  const __m128i swap_b32 = _mm_cvtepu16_epi32(swap_b);
+  // swap_b32 = [g f e 0] as i32
+  *x_sum_32 = _mm_add_epi32(*x_sum_32, sum_slli_a32);
+  *x_sum_32 = _mm_add_epi32(*x_sum_32, swap_b32);
+  // sum = [c+b+g a+f k+j+e i] as i32
+
+  // Also sum their squares
+  const __m128i slli_a_2 = _mm_madd_epi16(slli_a, slli_a);
+  const __m128i swap_b_2 = _mm_madd_epi16(swap_b, swap_b);
+  // slli_a_2 = [c2+b2 a2 k2+j2 i2]
+  // swap_b_2 = [0 0 g2+f2 e2]
+  const __m128i sum2 = _mm_hadd_epi32(slli_a_2, swap_b_2);
+  // sum2 = [0 g2+f2+e2 c2+b2+a2 k2+j2+i2]
+  *x2_sum_32 = _mm_add_epi32(*x2_sum_32, sum2);
+}
+
+void av1_get_horver_correlation_full_sse4_1(const int16_t *diff, int stride,
+                                            int width, int height, float *hcorr,
+                                            float *vcorr) {
+  // The following notation is used:
+  // x - current pixel
+  // y - right neighbour pixel
+  // z - below neighbour pixel
+  // w - down-right neighbour pixel
+  int64_t xy_sum = 0, xz_sum = 0;
+  int64_t x_sum = 0, x2_sum = 0;
+
+  // Process horizontal and vertical correlations through the body in 4x4
+  // blocks.  This excludes the final row and column and possibly one extra
+  // column depending how 3 divides into width and height
+  int32_t xy_tmp[4] = { 0 }, xz_tmp[4] = { 0 };
+  int32_t x_tmp[4] = { 0 }, x2_tmp[4] = { 0 };
+  __m128i xy_sum_32 = _mm_setzero_si128();
+  __m128i xz_sum_32 = _mm_setzero_si128();
+  __m128i x_sum_32 = _mm_setzero_si128();
+  __m128i x2_sum_32 = _mm_setzero_si128();
+  for (int i = 0; i <= height - 4; i += 3) {
+    for (int j = 0; j <= width - 4; j += 3) {
+      horver_correlation_4x4(&diff[i * stride + j], stride, &xy_sum_32,
+                             &xz_sum_32, &x_sum_32, &x2_sum_32);
+    }
+    xx_storeu_128(xy_tmp, xy_sum_32);
+    xx_storeu_128(xz_tmp, xz_sum_32);
+    xx_storeu_128(x_tmp, x_sum_32);
+    xx_storeu_128(x2_tmp, x2_sum_32);
+    xy_sum += (int64_t)xy_tmp[3] + xy_tmp[2] + xy_tmp[1];
+    xz_sum += (int64_t)xz_tmp[3] + xz_tmp[2] + xz_tmp[0];
+    x_sum += (int64_t)x_tmp[3] + x_tmp[2] + x_tmp[1] + x_tmp[0];
+    x2_sum += (int64_t)x2_tmp[2] + x2_tmp[1] + x2_tmp[0];
+    xy_sum_32 = _mm_setzero_si128();
+    xz_sum_32 = _mm_setzero_si128();
+    x_sum_32 = _mm_setzero_si128();
+    x2_sum_32 = _mm_setzero_si128();
+  }
+
+  // x_sum now covers every pixel except the final 1-2 rows and 1-2 cols
+  int64_t x_finalrow = 0, x_finalcol = 0, x2_finalrow = 0, x2_finalcol = 0;
+
+  // Do we have 2 rows remaining or just the one?  Note that width and height
+  // are powers of 2, so each modulo 3 must be 1 or 2.
+  if (height % 3 == 1) {  // Just horiz corrs on the final row
+    const int16_t x0 = diff[(height - 1) * stride];
+    x_sum += x0;
+    x_finalrow += x0;
+    x2_sum += x0 * x0;
+    x2_finalrow += x0 * x0;
+    for (int j = 0; j < width - 1; ++j) {
+      const int16_t x = diff[(height - 1) * stride + j];
+      const int16_t y = diff[(height - 1) * stride + j + 1];
+      xy_sum += x * y;
+      x_sum += y;
+      x2_sum += y * y;
+      x_finalrow += y;
+      x2_finalrow += y * y;
+    }
+  } else {  // Two rows remaining to do
+    const int16_t x0 = diff[(height - 2) * stride];
+    const int16_t z0 = diff[(height - 1) * stride];
+    x_sum += x0 + z0;
+    x2_sum += x0 * x0 + z0 * z0;
+    x_finalrow += z0;
+    x2_finalrow += z0 * z0;
+    for (int j = 0; j < width - 1; ++j) {
+      const int16_t x = diff[(height - 2) * stride + j];
+      const int16_t y = diff[(height - 2) * stride + j + 1];
+      const int16_t z = diff[(height - 1) * stride + j];
+      const int16_t w = diff[(height - 1) * stride + j + 1];
+
+      // Horizontal and vertical correlations for the penultimate row:
+      xy_sum += x * y;
+      xz_sum += x * z;
+
+      // Now just horizontal correlations for the final row:
+      xy_sum += z * w;
+
+      x_sum += y + w;
+      x2_sum += y * y + w * w;
+      x_finalrow += w;
+      x2_finalrow += w * w;
+    }
+  }
+
+  // Do we have 2 columns remaining or just the one?
+  if (width % 3 == 1) {  // Just vert corrs on the final col
+    const int16_t x0 = diff[width - 1];
+    x_sum += x0;
+    x_finalcol += x0;
+    x2_sum += x0 * x0;
+    x2_finalcol += x0 * x0;
+    for (int i = 0; i < height - 1; ++i) {
+      const int16_t x = diff[i * stride + width - 1];
+      const int16_t z = diff[(i + 1) * stride + width - 1];
+      xz_sum += x * z;
+      x_finalcol += z;
+      x2_finalcol += z * z;
+      // So the bottom-right elements don't get counted twice:
+      if (i < height - (height % 3 == 1 ? 2 : 3)) {
+        x_sum += z;
+        x2_sum += z * z;
+      }
+    }
+  } else {  // Two cols remaining
+    const int16_t x0 = diff[width - 2];
+    const int16_t y0 = diff[width - 1];
+    x_sum += x0 + y0;
+    x2_sum += x0 * x0 + y0 * y0;
+    x_finalcol += y0;
+    x2_finalcol += y0 * y0;
+    for (int i = 0; i < height - 1; ++i) {
+      const int16_t x = diff[i * stride + width - 2];
+      const int16_t y = diff[i * stride + width - 1];
+      const int16_t z = diff[(i + 1) * stride + width - 2];
+      const int16_t w = diff[(i + 1) * stride + width - 1];
+
+      // Horizontal and vertical correlations for the penultimate col:
+      // Skip these on the last iteration of this loop if we also had two
+      // rows remaining, otherwise the final horizontal and vertical correlation
+      // get erroneously processed twice
+      if (i < height - 2 || height % 3 == 1) {
+        xy_sum += x * y;
+        xz_sum += x * z;
+      }
+
+      x_finalcol += w;
+      x2_finalcol += w * w;
+      // So the bottom-right elements don't get counted twice:
+      if (i < height - (height % 3 == 1 ? 2 : 3)) {
+        x_sum += z + w;
+        x2_sum += z * z + w * w;
+      }
+
+      // Now just vertical correlations for the final column:
+      xz_sum += y * w;
+    }
+  }
+
+  // Calculate the simple sums and squared-sums
+  int64_t x_firstrow = 0, x_firstcol = 0;
+  int64_t x2_firstrow = 0, x2_firstcol = 0;
+
+  for (int j = 0; j < width; ++j) {
+    x_firstrow += diff[j];
+    x2_firstrow += diff[j] * diff[j];
+  }
+  for (int i = 0; i < height; ++i) {
+    x_firstcol += diff[i * stride];
+    x2_firstcol += diff[i * stride] * diff[i * stride];
+  }
+
+  int64_t xhor_sum = x_sum - x_finalcol;
+  int64_t xver_sum = x_sum - x_finalrow;
+  int64_t y_sum = x_sum - x_firstcol;
+  int64_t z_sum = x_sum - x_firstrow;
+  int64_t x2hor_sum = x2_sum - x2_finalcol;
+  int64_t x2ver_sum = x2_sum - x2_finalrow;
+  int64_t y2_sum = x2_sum - x2_firstcol;
+  int64_t z2_sum = x2_sum - x2_firstrow;
+
+  const float num_hor = (float)(height * (width - 1));
+  const float num_ver = (float)((height - 1) * width);
+
+  const float xhor_var_n = x2hor_sum - (xhor_sum * xhor_sum) / num_hor;
+  const float xver_var_n = x2ver_sum - (xver_sum * xver_sum) / num_ver;
+
+  const float y_var_n = y2_sum - (y_sum * y_sum) / num_hor;
+  const float z_var_n = z2_sum - (z_sum * z_sum) / num_ver;
+
+  const float xy_var_n = xy_sum - (xhor_sum * y_sum) / num_hor;
+  const float xz_var_n = xz_sum - (xver_sum * z_sum) / num_ver;
+
+  if (xhor_var_n > 0 && y_var_n > 0) {
+    *hcorr = xy_var_n / sqrtf(xhor_var_n * y_var_n);
+    *hcorr = *hcorr < 0 ? 0 : *hcorr;
+  } else {
+    *hcorr = 1.0;
+  }
+  if (xver_var_n > 0 && z_var_n > 0) {
+    *vcorr = xz_var_n / sqrtf(xver_var_n * z_var_n);
+    *vcorr = *vcorr < 0 ? 0 : *vcorr;
+  } else {
+    *vcorr = 1.0;
+  }
+}
diff --git a/test/horver_correlation_test.cc b/test/horver_correlation_test.cc
new file mode 100644
index 0000000..5fca4b6
--- /dev/null
+++ b/test/horver_correlation_test.cc
@@ -0,0 +1,149 @@
+/*
+ * Copyright (c) 2018, Alliance for Open Media. All rights reserved
+ *
+ * This source code is subject to the terms of the BSD 2 Clause License and
+ * the Alliance for Open Media Patent License 1.0. If the BSD 2 Clause License
+ * was not distributed with this source code in the LICENSE file, you can
+ * obtain it at www.aomedia.org/license/software. If the Alliance for Open
+ * Media Patent License 1.0 was not distributed with this source code in the
+ * PATENTS file, you can obtain it at www.aomedia.org/license/patent.
+ */
+
+#include <vector>
+
+#include "third_party/googletest/src/googletest/include/gtest/gtest.h"
+
+#include "test/function_equivalence_test.h"
+#include "test/register_state_check.h"
+
+#include "config/aom_config.h"
+#include "config/aom_dsp_rtcd.h"
+#include "config/av1_rtcd.h"
+
+#include "aom/aom_integer.h"
+#include "av1/encoder/rdopt.h"
+
+using libaom_test::FunctionEquivalenceTest;
+
+namespace {
+typedef void (*horver_Func)(const int16_t *diff, int stride, int w, int h,
+                            float *hcorr, float *vcorr);
+
+typedef libaom_test::FuncParam<horver_Func> TestFuncs;
+
+typedef ::testing::tuple<const horver_Func> HorverTestParam;
+
+class HorverTest : public ::testing::TestWithParam<HorverTestParam> {
+ public:
+  virtual void SetUp() {
+    data_buf = (int16_t *)aom_malloc(MAX_SB_SQUARE * sizeof(int16_t));
+    target_func_ = GET_PARAM(0);
+  }
+  virtual void TearDown() { aom_free(data_buf); }
+  void runHorverTest(void);
+  void runHorverTest_ExtremeValues(void);
+  void runHorverSpeedTest(int run_times);
+
+ private:
+  horver_Func target_func_;
+  ACMRandom rng_;
+  int16_t *data_buf;
+};
+
+void HorverTest::runHorverTest(void) {
+  for (int block_size = 0; block_size < BLOCK_SIZES_ALL; block_size++) {
+    const int w = block_size_wide[block_size];
+    const int h = block_size_high[block_size];
+    for (int iter = 0; iter < 1000 && !HasFatalFailure(); ++iter) {
+      float hcorr_ref = 0.0, vcorr_ref = 0.0;
+      float hcorr_test = 0.0, vcorr_test = 0.0;
+
+      for (int i = 0; i < MAX_SB_SQUARE; ++i) {
+        data_buf[i] = rng_.Rand16() % (1 << 12);
+      }
+
+      av1_get_horver_correlation_full_c(data_buf, MAX_SB_SIZE, w, h, &hcorr_ref,
+                                        &vcorr_ref);
+
+      target_func_(data_buf, MAX_SB_SIZE, w, h, &hcorr_test, &vcorr_test);
+
+      ASSERT_LE(fabs(hcorr_ref - hcorr_test), 1e-6)
+          << "hcorr incorrect (" << w << "x" << h << ")";
+      ASSERT_LE(fabs(vcorr_ref - vcorr_test), 1e-6)
+          << "vcorr incorrect (" << w << "x" << h << ")";
+    }
+    //    printf("(%3dx%-3d) passed\n", w, h);
+  }
+}
+
+void HorverTest::runHorverSpeedTest(int run_times) {
+  for (int i = 0; i < MAX_SB_SQUARE; ++i) {
+    data_buf[i] = rng_.Rand16() % (1 << 12);
+  }
+
+  for (int block_size = 0; block_size < BLOCK_SIZES_ALL; block_size++) {
+    const int w = block_size_wide[block_size];
+    const int h = block_size_high[block_size];
+    float hcorr_ref = 0.0, vcorr_ref = 0.0;
+    float hcorr_test = 0.0, vcorr_test = 0.0;
+
+    aom_usec_timer timer;
+    aom_usec_timer_start(&timer);
+    for (int i = 0; i < run_times; ++i) {
+      av1_get_horver_correlation_full_c(data_buf, MAX_SB_SIZE, w, h, &hcorr_ref,
+                                        &vcorr_ref);
+    }
+    aom_usec_timer_mark(&timer);
+    const double time1 = static_cast<double>(aom_usec_timer_elapsed(&timer));
+    aom_usec_timer_start(&timer);
+    for (int i = 0; i < run_times; ++i) {
+      target_func_(data_buf, MAX_SB_SIZE, w, h, &hcorr_test, &vcorr_test);
+    }
+    aom_usec_timer_mark(&timer);
+    const double time2 = static_cast<double>(aom_usec_timer_elapsed(&timer));
+
+    printf("%3dx%-3d:%7.2f/%7.2fns (%3.2f)\n", w, h, time1, time2,
+           time1 / time2);
+  }
+}
+
+void HorverTest::runHorverTest_ExtremeValues(void) {
+  for (int i = 0; i < MAX_SB_SQUARE; ++i) {
+    // Most of get_horver_test is squaring and summing, so simply saturating
+    // the whole buffer is mostly likely to cause an overflow.
+    data_buf[i] = (1 << 12) - 1;
+  }
+
+  for (int block_size = 0; block_size < BLOCK_SIZES_ALL; block_size++) {
+    const int w = block_size_wide[block_size];
+    const int h = block_size_high[block_size];
+    float hcorr_ref = 0.0, vcorr_ref = 0.0;
+    float hcorr_test = 0.0, vcorr_test = 0.0;
+
+    av1_get_horver_correlation_full_c(data_buf, MAX_SB_SIZE, w, h, &hcorr_ref,
+                                      &vcorr_ref);
+    target_func_(data_buf, MAX_SB_SIZE, w, h, &hcorr_test, &vcorr_test);
+
+    ASSERT_LE(fabs(hcorr_ref - hcorr_test), 1e-6) << "hcorr incorrect";
+    ASSERT_LE(fabs(vcorr_ref - vcorr_test), 1e-6) << "vcorr incorrect";
+  }
+}
+
+TEST_P(HorverTest, RandomValues) { runHorverTest(); }
+
+TEST_P(HorverTest, ExtremeValues) { runHorverTest_ExtremeValues(); }
+
+TEST_P(HorverTest, DISABLED_Speed) { runHorverSpeedTest(100000); }
+
+#if HAVE_SSE4_1
+INSTANTIATE_TEST_CASE_P(
+    SSE4_1, HorverTest,
+    ::testing::Values(av1_get_horver_correlation_full_sse4_1));
+#endif  // HAVE_SSE4_1
+
+#if HAVE_AVX2
+INSTANTIATE_TEST_CASE_P(
+    AVX2, HorverTest, ::testing::Values(av1_get_horver_correlation_full_avx2));
+#endif  // HAVE_AVX2
+
+}  // namespace
diff --git a/test/test.cmake b/test/test.cmake
index 513a0f1..d15e580 100644
--- a/test/test.cmake
+++ b/test/test.cmake
@@ -184,6 +184,7 @@
               "${AOM_ROOT}/test/error_block_test.cc"
               "${AOM_ROOT}/test/fft_test.cc"
               "${AOM_ROOT}/test/fwht4x4_test.cc"
+              "${AOM_ROOT}/test/horver_correlation_test.cc"
               "${AOM_ROOT}/test/masked_sad_test.cc"
               "${AOM_ROOT}/test/masked_variance_test.cc"
               "${AOM_ROOT}/test/motion_vector_test.cc"