diff --git a/python/lsst/obs/subaru/isr.py b/python/lsst/obs/subaru/isr.py
index ce7449dd..51f3b3bd 100644
--- a/python/lsst/obs/subaru/isr.py
+++ b/python/lsst/obs/subaru/isr.py
@@ -44,6 +44,9 @@ import lsst.afw.display.ds9 as ds9
 from lsst.obs.hsc.vignette import VignetteConfig
 from lsst.afw.geom.polygon import Polygon
 
+from .straylightRemoval import waveletCompression
+from .straylightRemoval import SCAGFunctions
+
 
 class QaFlatnessConfig(pexConfig.Config):
     meshX = pexConfig.Field(
@@ -108,6 +111,16 @@ class SubaruIsrConfig(IsrTask.ConfigClass):
         doc="Correct for crosstalk",
         default=True,
     )
+    doSubtractStraylight = pexConfig.Field(
+        dtype=bool,
+        doc="subtract straylights?",
+        default=True,
+    )
+    doKawanomotoCorrection = pexConfig.Field(
+        dtype=bool,
+        doc="do Kawanomoto correction in straylights subtraction?",
+        default=True,
+    )
     doApplyGains = pexConfig.Field(
         dtype=bool,
         doc="""Correct the amplifiers for their gains
@@ -305,6 +318,10 @@ class SubaruIsrTask(IsrTask):
         if self.config.doDark:
             darkExposure = self.getIsrExposure(sensorRef, "dark")
             self.darkCorrection(ccdExposure, darkExposure)
+
+        if self.config.doSubtractStraylight:
+            self.subtractStraylight(sensorRef, ccdExposure, self.config)
+
         if self.config.doFlat:
             flatExposure = self.getIsrExposure(sensorRef, "flat")
             self.flatCorrection(ccdExposure, flatExposure)
@@ -587,6 +604,45 @@ class SubaruIsrTask(IsrTask):
         self.log.info("Setting rough magnitude zero point: %f" % (2.5*math.log10(fluxMag0*expTime),))
         exposure.getCalib().setFluxMag0(fluxMag0*expTime)
 
+    def subtractStraylight(self, sensorRef, exposure, config):
+        filterName = exposure.getFilter().getName()
+        if filterName != 'y':
+            return
+
+        # The LEDs that are causing the Y straylight have not been covered yet (on 2017-11-27),
+        # but they will be covered in near future.
+        # Once the LEDs are covered, we will have to uncomment the following statement:
+        #
+        # if (ccdExposure is newer than a certain date):
+        #     return
+
+        ybackground = waveletCompression.ybackground
+        ybackground.config.data_dir = os.path.join(os.path.dirname(__file__), 'straylightRemoval', 'dictionary', filterName)
+
+        dictionaryEnv = 'STRAYLIGHT_DICTIONARY'
+        if dictionaryEnv in os.environ:
+            ybackground.config.data_dir = os.environ[dictionaryEnv]
+
+        header = sensorRef.get('raw_md')
+        ccd = sensorRef.dataId['ccd']
+
+        if config.doKawanomotoCorrection:
+            angleStart, angleEnd = SCAGFunctions.inrStartEnd(header)
+
+            self.log.info(
+                "(INR-STR, INR-END) = ({:g}, {:g}) (FITS header says ({:g}, {:g}))."
+                .format(angleStart, angleEnd, header.getDouble('INR-STR'), header.getDouble('INR-END'))
+            )
+
+        else:
+            angleStart = angleEnd = header.getDouble('INR-STR')
+
+        expTime = exposure.getInfo().getVisitInfo().getExposureTime()
+        straylight = ybackground.get_ybackground(ccd, angleStart, None if angleStart == angleEnd else angleEnd)
+        array = exposure.getMaskedImage().getImage().getArray()
+        array -= expTime * straylight
+        return exposure
+
 
 class SuprimecamIsrConfig(SubaruIsrConfig):
 
diff --git a/python/lsst/obs/subaru/straylightRemoval/SCAGFunctions.py b/python/lsst/obs/subaru/straylightRemoval/SCAGFunctions.py
new file mode 100644
index 00000000..4352f28e
--- /dev/null
+++ b/python/lsst/obs/subaru/straylightRemoval/SCAGFunctions.py
@@ -0,0 +1,234 @@
+#!/usr/bin/env python
+
+# Copyright (C) 2017  HSC Software Team
+# Copyright (C) 2017  Satoshi Kawanomoto
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+import astropy.io.fits as pyfits
+import numpy as np
+import sys
+
+### fixed parameters
+ltt_d     =   19.82556 # dome latitude in degree
+lng_d     = -155.47611 # dome longitude in degree
+mjd_J2000 =    51544.5 # mjd at J2000.0 (2000/01/01.5)
+
+# refraction index of air
+# T=273.15[K], P=600[hPa], Pw=1.5[hPa], lambda=0.55[um]
+air_idx = 1.0 + 1.7347e-04
+
+# scale height of air
+air_sh = 0.00130
+
+def _mjd2jc2000(mjd):
+    """convert mjd to Julian century (J2000.0 origin)"""
+    jc2000 = (mjd - mjd_J2000) / 36525.0
+    return jc2000
+
+def _precessionMatrix(jc2000):
+    """create precession matrix at the given time in Julian century"""
+    zeta_A  = np.deg2rad((2306.2181*jc2000 + 0.30188*jc2000**2.0 + 0.017998*jc2000**3.0)/3600.0)
+    z_A     = np.deg2rad((2306.2181*jc2000 + 1.09468*jc2000**2.0 + 0.018203*jc2000**3.0)/3600.0)
+    theta_A = np.deg2rad((2004.3109*jc2000 - 0.42665*jc2000**2.0 - 0.041833*jc2000**3.0)/3600.0)
+    precMat = np.matrix([[+np.cos(zeta_A)*np.cos(theta_A)*np.cos(z_A)-np.sin(zeta_A)*np.sin(z_A), \
+                          -np.sin(zeta_A)*np.cos(theta_A)*np.cos(z_A)-np.cos(zeta_A)*np.sin(z_A), \
+                          -np.sin(theta_A)*np.cos(z_A)],
+                         [+np.cos(zeta_A)*np.cos(theta_A)*np.sin(z_A)+np.sin(zeta_A)*np.cos(z_A), \
+                          -np.sin(zeta_A)*np.cos(theta_A)*np.sin(z_A)+np.cos(zeta_A)*np.cos(z_A), \
+                          -np.sin(theta_A)*np.sin(z_A)],
+                         [+np.cos(zeta_A)*np.sin(theta_A)                                       , \
+                          -np.sin(zeta_A)*np.sin(theta_A)                                       , \
+                          +np.cos(theta_A)]])
+    return precMat
+
+def _mjd2gmst(mjd):
+    """convert mjd to GMST(Greenwich mean sidereal time)"""
+    mjd_f = mjd%1
+    jc2000 = _mjd2jc2000(mjd)
+    gmst_s = (6.0*3600.0+41.0*60.0+50.54841)\
+              + 8640184.812866*jc2000 + 0.093104*jc2000**2.0 - 0.0000062*jc2000**3.0\
+              + mjd_f*86400.0
+    gmst_d = (gmst_s%86400)/240.0
+    return gmst_d
+
+def _gmst2lmst(gmst_d):
+    """convert GMST to LMST(mean local sidereal time)"""
+    lmst_d = (gmst_d + lng_d)%360
+    return lmst_d
+
+def _sph2vec(ra_d, de_d):
+    """convert spherical coordinate to the Cartesian coordinates (vector)"""
+    ra_r = np.deg2rad(ra_d)
+    de_r = np.deg2rad(de_d)
+    vec  = np.array([[np.cos(ra_r)*np.cos(de_r)],
+                     [np.sin(ra_r)*np.cos(de_r)],
+                     [             np.sin(de_r)]])
+    return vec
+
+def _vec2sph(vec):
+    """convert the Cartesian coordinates vector to shperical coordinates"""
+    ra_r = np.arctan2(vec[1, 0], vec[0, 0])
+    de_r = np.arcsin(vec[2, 0])
+    ra_d = np.rad2deg(ra_r)
+    de_d = np.rad2deg(de_r)
+    return ra_d, de_d
+
+def _ra2ha(ra_d, lst_d):
+    """convert right ascension to hour angle at given LST"""
+    ha_d = (lst_d - ra_d)%360
+    return ha_d
+
+def _eq2hz(ha_d, de_d):
+    """convert equatorial coordinates to the horizontal coordinates"""
+    ltt_r = np.deg2rad(ltt_d)
+    ha_r  = np.deg2rad(ha_d)
+    de_r  = np.deg2rad(de_d)
+    zd_r  = np.arccos(+np.sin(ltt_r)*np.sin(de_r)+np.cos(ltt_r)*np.cos(de_r)*np.cos(ha_r))
+    az_r  = np.arctan2(+np.cos(de_r)*np.sin(ha_r),
+                       -np.cos(ltt_r)*np.sin(de_r)+np.sin(ltt_r)*np.cos(de_r)*np.cos(ha_r))
+    zd_d  = np.rad2deg(zd_r)
+    az_d  = np.rad2deg(az_r)
+    al_d  = 90.0 - zd_d
+    return al_d, az_d
+
+def _air_idx():
+    """return the air refraction index"""
+    return air_idx
+
+def _atm_ref(al_d):
+    """return the atmospheric refraction at given altitude"""
+    if al_d > 20.0:
+        zd_r = np.deg2rad(90.0-al_d)
+    else:
+        zd_r = np.deg2rad(70.0)
+    r0 = _air_idx()-1.0
+    sh = air_sh
+    R0 = (1.0-sh)*r0-sh*r0**2/2.0+sh**2*r0*2.0
+    R1 = r0**2/2.0+r0**3/6.0-sh*r0-sh*r0**2*11.0/4.0+sh**2*r0*5.0
+    R2 = r0**3-sh*r0**2*9.0/4.0+sh**2*r0*3.0
+    R = R0*np.tan(zd_r)+R1*(np.tan(zd_r))**3+R2*(np.tan(zd_r))**5
+    return np.rad2deg(R)
+
+def _mal2aal(mal_d):
+    """convert mean altitude to apparent altitude"""
+    aal_d = mal_d + _atm_ref(mal_d)
+    return aal_d
+
+def _pos2adt(al_t_d, al_s_d, delta_az_d):
+    """convert altitudes of telescope and star and relative azimuth to angular distance and position angle"""
+    zd_t_r = np.deg2rad(90.0-al_t_d)
+    zd_s_r = np.deg2rad(90.0-al_s_d)
+    daz_r  = np.deg2rad(delta_az_d)
+
+    ad_r   = np.arccos(np.cos(zd_t_r)*np.cos(zd_s_r)+np.sin(zd_t_r)*np.sin(zd_s_r)*np.cos(daz_r))
+
+    if ad_r > 0.0:
+        pa_r = np.arcsin(np.sin(zd_s_r)*np.sin(daz_r)/np.sin(ad_r))
+    else:
+        pa_r = 0.0
+    ad_d = np.rad2deg(ad_r)
+    pa_d = np.rad2deg(pa_r)
+    if (zd_t_r < zd_s_r):
+        pa_d = 180.0 - pa_d
+
+    return ad_d, pa_d
+
+def _addpad2xy(ang_dist_d, p_ang_d, inr_d):
+    """convert angular distance, position angle, and instrument rotator angle to position on the cold plate"""
+    t = 90.0-(p_ang_d-inr_d)
+    x = np.cos(np.deg2rad(t))
+    y = np.sin(np.deg2rad(t))
+    return x, y
+
+def _gsCPposNorth(ra_t_d, de_t_d, inr_d, mjd):
+    jc2000 = _mjd2jc2000(mjd)
+    pm     = _precessionMatrix(jc2000)
+
+    vt      = _sph2vec(ra_t_d, de_t_d)
+    vt_mean = np.dot(pm, vt)
+
+    (mean_ra_t_d, mean_de_t_d) = _vec2sph(vt_mean)
+    mean_ra_s_d = mean_ra_t_d
+    mean_de_s_d = mean_de_t_d+0.75
+
+    gmst_d = _mjd2gmst(mjd)
+    lmst_d = _gmst2lmst(gmst_d)
+
+    mean_ha_t_d = _ra2ha(mean_ra_t_d, lmst_d)
+    mean_ha_s_d = _ra2ha(mean_ra_s_d, lmst_d)
+
+    (mean_al_t_d, mean_az_t_d) =  _eq2hz(mean_ha_t_d, mean_de_t_d)
+    (mean_al_s_d, mean_az_s_d) =  _eq2hz(mean_ha_s_d, mean_de_s_d)
+
+    apparent_al_t_d = _mal2aal(mean_al_t_d)
+    apparent_al_s_d = _mal2aal(mean_al_s_d)
+
+    delta_az_d = mean_az_s_d - mean_az_t_d
+
+    (ang_dist_d, p_ang_d) = _pos2adt(apparent_al_t_d, apparent_al_s_d, delta_az_d)
+
+    (x, y) = _addpad2xy(ang_dist_d, p_ang_d, inr_d)
+    return x, y
+
+def _getDataArrayFromFITSFileWithHeader(pathToFITSFile):
+    """return array of pixel data"""
+    fitsfile   = pyfits.open(pathToFITSFile)
+    dataArray  = fitsfile[0].data
+    fitsHeader = fitsfile[0].header
+    fitsfile.close()
+    return dataArray, fitsHeader
+
+def _minorArc(angle1, angle2):
+    """e.g. input (-179, 179) -> output (-179, -181)"""
+
+    angle1 = (angle1 + 180.0) % 360 - 180.0
+    angle2 = (angle2 + 180.0) % 360 - 180.0
+
+    if angle1 < angle2:
+        if angle2 - angle1 > 180.0:
+            angle2 -= 360.0
+    elif angle2 < angle1:
+        if angle1 - angle2 > 180.0:
+            angle1 -= 360.0
+
+    # Try to place [angle1, angle2] within [-270, +270]
+
+    if min(angle1, angle2) < -270.0:
+        angle1 += 360.0
+        angle2 += 360.0
+    if max(angle1, angle2) > 270.0:
+        angle1 -= 360.0
+        angle2 -= 360.0
+
+    return angle1, angle2
+
+def inrStartEnd(header):
+    # header: PropertyList http://doxygen.lsst.codes/stack/doxygen/x_masterDoxyDoc/classlsst_1_1daf_1_1base_1_1_property_list.html
+
+    inst_pa = header.getDouble('INST-PA')
+    ra_t_d = header.getDouble('CRVAL1')
+    de_t_d = header.getDouble('CRVAL2')
+
+    mjd_str = header.getDouble('MJD-STR')
+    mjd_end = header.getDouble('MJD-END')
+
+    inr_d = 0.00
+
+    (x, y) = _gsCPposNorth(ra_t_d, de_t_d, inr_d, mjd_str)
+    x_inr_str = 90.0-np.rad2deg(np.arctan2(y, x))+inst_pa
+    (x, y) = _gsCPposNorth(ra_t_d, de_t_d, inr_d, mjd_end)
+    x_inr_end =  90.0-np.rad2deg(np.arctan2(y, x))+inst_pa
+
+    return _minorArc(x_inr_str, x_inr_end)
diff --git a/python/lsst/obs/subaru/straylightRemoval/waveletCompression/__init__.py b/python/lsst/obs/subaru/straylightRemoval/waveletCompression/__init__.py
new file mode 100644
index 00000000..69ab2b3e
--- /dev/null
+++ b/python/lsst/obs/subaru/straylightRemoval/waveletCompression/__init__.py
@@ -0,0 +1 @@
+from .ybackground import *
diff --git a/python/lsst/obs/subaru/straylightRemoval/waveletCompression/config.py b/python/lsst/obs/subaru/straylightRemoval/waveletCompression/config.py
new file mode 100644
index 00000000..1bf9f64e
--- /dev/null
+++ b/python/lsst/obs/subaru/straylightRemoval/waveletCompression/config.py
@@ -0,0 +1,4 @@
+import os.path
+
+rawdata_dir = None
+data_dir = os.path.join(os.path.dirname(__file__), '..', 'dictionary', 'y')
diff --git a/python/lsst/obs/subaru/straylightRemoval/waveletCompression/dwt.py b/python/lsst/obs/subaru/straylightRemoval/waveletCompression/dwt.py
new file mode 100644
index 00000000..64f1e875
--- /dev/null
+++ b/python/lsst/obs/subaru/straylightRemoval/waveletCompression/dwt.py
@@ -0,0 +1,421 @@
+# ybackground
+# Copyright (C) 2017  Sogo Mineo
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+#----------------------------------------------------------------------------
+# CDF 9/7 wavelet transformation
+# See
+#   T.800: JPEG2000 image coding system: Core coding system
+#       https://www.itu.int/rec/T-REC-T.800
+#   OpenJPEG
+#       https://github.com/uclouvain/openjpeg
+
+from __future__ import division, absolute_import, print_function
+
+import numpy
+
+import itertools
+
+# From Table F.4 in ITU-T Rec. T.800
+# Definition of lifting parameters for the 9-7 irreversible filter
+alpha = -1.586134342059924
+beta  = -0.052980118572961
+gamma =  0.882911075530934
+delta =  0.443506852043971
+# K     =  1.230174104914001
+
+# "K" determines the normalization of biorthogonal filters.
+# JPEG2000 uses the normalization:
+#    H(0) = \sum_i h[i] = 1, \tilde{H}(0) = \sum_i \tilde{h}[i] = 2
+# where h is the analysis low-pass filter, and \tilde{h} is the synthesis low-pass filter.
+# The value K = 1.230174104914001 is for this normalization.
+#
+# We use the normalization:
+#    H(0) = \sqrt{2}, \tilde{H}(0) = \sqrt{2}.
+# For this normalization, K must be the following value.
+K     =  0.86986445162477855
+
+def cdf_9_7(data, level):
+    """
+    Multi-dimensional forward wavelet transformation
+    with the Cohen-Daubechies-Feauveau wavelet 9/7
+
+    @param data (numpy.array)
+        N-dimensional array (N >= 1). This array will be *destroyed*.
+        If you want `data` kept, copy it before calling this function:
+            cdf_9_7(numpy.copy(data), level)
+
+    @param level (int or tuple)
+        Level of wavelet analysis.
+        If a tuple is given, its length must agree with the rank of `data`,
+        and level[i] is to be the level for the i-th axis.
+
+    @return (numpy.array)
+        Result of multiresolution analysis.
+    """
+    data = numpy.asarray(data, dtype=float)
+    if numpy.asarray(level, dtype=int).shape == ():
+        level = itertools.repeat(level)
+
+    for axis, lev in zip(range(len(data.shape)), level):
+        data = cdf_9_7_1d(data, lev, axis)
+
+    return data
+
+def cdf_9_7_1d(data, level, axis=0):
+    """
+    One-dimensional forward wavelet transformation.
+
+    @param data (numpy.array)
+        N-dimensional array (N >= 1). This array will be *destroyed*.
+
+    @param level (int)
+        Level of wavelet analysis.
+
+    @param axis (int)
+
+    @return (numpy.array)
+        Result of multiresolution analysis.
+    """
+    data = numpy.asarray(data, dtype=float)
+    data = numpy.swapaxes(data, 0, axis)
+    size = len(data)
+
+    for i in range(level):
+        data[:size] = _cdf_9_7_1d_level1(data[:size])
+        size = (size + 1) // 2
+
+    return numpy.swapaxes(data, 0, axis)
+
+def _cdf_9_7_1d_level1(data):
+    """
+    Level 1 one-dimensional forward wavelet transformation.
+    The 0th axis is used in the transformation.
+
+    @param data (numpy.array)
+        N-dimensional array (N >= 1). This array will be *destroyed*.
+
+    @return (numpy.array)
+        Result of transformation.
+    """
+    size = len(data)
+    if size <= 1:
+        return data
+
+    iseven = (size % 2 == 0)
+
+    # From equations (F-11) in ITU-T Rec. T.800
+
+    # predict (1)
+    data[1:-1:2] += alpha * (data[0:-2:2] + data[2::2])
+    if iseven:
+        data[-1] += 2*alpha * data[-2]
+
+    # update (1)
+    data[0] += 2*beta * data[1]
+    data[2:-1:2] += beta * (data[1:-2:2] + data[3::2])
+    if not iseven:
+        data[-1] += 2*beta * data[-2]
+
+    # predict (2)
+    data[1:-1:2] += gamma * (data[0:-2:2] + data[2::2])
+    if iseven:
+        data[-1] += 2*gamma * data[-2]
+
+    # update (2)
+    data[0] += 2*delta * data[1]
+    data[2:-1:2] += delta * (data[1:-2:2] + data[3::2])
+    if not iseven:
+        data[-1] += 2*delta * data[-2]
+
+    # de-interleave
+    scaling_size = (size + 1) // 2
+    ret = numpy.empty_like(data)
+    ret[:scaling_size] = (1/K) * data[0::2]
+    ret[scaling_size:] = K     * data[1::2]
+
+    return ret
+
+
+def icdf_9_7(data, level):
+    """
+    Multi-dimensional backword wavelet transformation
+    with the Cohen-Daubechies-Feauveau wavelet 9/7
+
+    @param data (numpy.array)
+        N-dimensional array (N >= 1). This array will be *destroyed*.
+        If you want `data` kept, copy it before calling this function:
+            icdf_9_7(numpy.copy(data), level)
+
+    @param level (int or tuple)
+        Level of wavelet analysis.
+        If a tuple is given, its length must agree with the rank of `data`,
+        and level[i] is to be the level for the i-th axis.
+
+    @return (numpy.array)
+        Result of multiresolution synthesis.
+    """
+    data = numpy.asarray(data, dtype=float)
+    if numpy.asarray(level, dtype=int).shape == ():
+        level = itertools.repeat(level)
+
+    for axis, lev in zip(range(len(data.shape)), level):
+        data = icdf_9_7_1d(data, lev, axis)
+
+    return data
+
+def icdf_9_7_1d(data, level, axis=0):
+    """
+    One-dimensional backword wavelet transformation.
+    The 0th axis is used in the transformation.
+
+    @param data (numpy.array)
+        N-dimensional array (N >= 1). This array will be *destroyed*.
+
+    @param level (int)
+        Level of wavelet analysis.
+
+    @param axis (int)
+
+    @return (numpy.array)
+        Result of multiresolution synthesis.
+    """
+    data = numpy.asarray(data, dtype=float)
+    data = numpy.swapaxes(data, 0, axis)
+    size = len(data)
+    sizes = []
+    for i in range(level):
+        sizes.append(size)
+        size = (size + 1) // 2
+
+    for size in reversed(sizes):
+        data[:size] = _icdf_9_7_1d_level1(data[:size])
+
+    return numpy.swapaxes(data, 0, axis)
+
+def _icdf_9_7_1d_level1(data):
+    """
+    Level 1 one-dimensional backword wavelet transformation.
+    The 0th axis is used in the transformation.
+
+    @param data (numpy.array)
+        N-dimensional array (N >= 1). This array will be *destroyed*.
+
+    @return (numpy.array)
+        Result of transformation.
+    """
+    size = len(data)
+    if size <= 1:
+        return data
+
+    iseven = (size % 2 == 0)
+
+    # interleave
+    scaling_size = (size + 1) // 2
+    ret = data
+    data = numpy.empty_like(ret)
+    data[0::2] = K     * ret[:scaling_size]
+    data[1::2] = (1/K) * ret[scaling_size:]
+
+    # From equations (F-7) in ITU-T Rec. T.800
+
+    # update (2)
+    data[0] -= 2*delta * data[1]
+    data[2:-1:2] -= delta * (data[1:-2:2] + data[3::2])
+    if not iseven:
+        data[-1] -= 2*delta * data[-2]
+
+    # predict (2)
+    data[1:-1:2] -= gamma * (data[0:-2:2] + data[2::2])
+    if iseven:
+        data[-1] -= 2*gamma * data[-2]
+
+    # update (1)
+    data[0] -= 2*beta * data[1]
+    data[2:-1:2] -= beta * (data[1:-2:2] + data[3::2])
+    if not iseven:
+        data[-1] -= 2*beta * data[-2]
+
+    # predict (1)
+    data[1:-1:2] -= alpha * (data[0:-2:2] + data[2::2])
+    if iseven:
+        data[-1] -= 2*alpha * data[-2]
+
+    return data
+
+
+def periodic_cdf_9_7_1d(data, level, axis=0):
+    """
+    One-dimensional forward wavelet transformation.
+
+    @param data (numpy.array)
+        N-dimensional array (N >= 1). This array will be *destroyed*.
+        The size of this data must be divisible by 2**level.
+
+    @param level (int)
+        Level of wavelet analysis.
+
+    @param axis (int)
+
+    @return (numpy.array)
+        Result of multiresolution analysis.
+    """
+    data = numpy.asarray(data, dtype=float)
+    data = numpy.swapaxes(data, 0, axis)
+    size = len(data)
+
+    if size % (2**level) != 0:
+        raise ValueError("Size must be divisible by 2**level.")
+
+    for i in range(level):
+        data[:size] = _periodic_cdf_9_7_1d_level1(data[:size])
+        size = (size + 1) // 2
+
+    return numpy.swapaxes(data, 0, axis)
+
+def _periodic_cdf_9_7_1d_level1(data):
+    """
+    Level 1 one-dimensional forward wavelet transformation.
+    The 0th axis is used in the transformation.
+
+    @param data (numpy.array)
+        N-dimensional array (N >= 1). This array will be *destroyed*.
+
+    @return (numpy.array)
+        Result of transformation.
+    """
+    size = len(data)
+    if size <= 1:
+        return data
+
+    assert(size % 2 == 0)
+
+    # From equations (F-11) in ITU-T Rec. T.800
+
+    # predict (1)
+    data[1:-1:2] += alpha * (data[0:-2:2] + data[2::2])
+    data[-1] += alpha * (data[0] + data[-2])
+
+    # update (1)
+    data[0] += beta * (data[1] + data[-1])
+    data[2:-1:2] += beta * (data[1:-2:2] + data[3::2])
+
+    # predict (2)
+    data[1:-1:2] += gamma * (data[0:-2:2] + data[2::2])
+    data[-1] += gamma * (data[0] + data[-2])
+
+    # update (2)
+    data[0] += delta * (data[1] + data[-1])
+    data[2:-1:2] += delta * (data[1:-2:2] + data[3::2])
+
+    # de-interleave
+    scaling_size = (size + 1) // 2
+    ret = numpy.empty_like(data)
+    ret[:scaling_size] = (1/K) * data[0::2]
+    ret[scaling_size:] = K     * data[1::2]
+
+    return ret
+
+
+def periodic_icdf_9_7_1d(data, level, axis=0):
+    """
+    One-dimensional backword wavelet transformation.
+
+    @param data (numpy.array)
+        N-dimensional array (N >= 1). This array will be *destroyed*.
+        The size of this data must be divisible by 2**level.
+
+    @param level (int)
+        Level of wavelet analysis.
+
+    @param axis (int)
+
+    @return (numpy.array)
+        Result of multiresolution synthesis.
+    """
+    data = numpy.asarray(data, dtype=float)
+    data = numpy.swapaxes(data, 0, axis)
+    size = len(data)
+
+    if size % (2**level) != 0:
+        raise ValueError("Size must be divisible by 2**level.")
+
+    sizes = []
+    for i in range(level):
+        sizes.append(size)
+        size = (size + 1) // 2
+
+    for size in reversed(sizes):
+        data[:size] = _periodic_icdf_9_7_1d_level1(data[:size])
+
+    return numpy.swapaxes(data, 0, axis)
+
+def _periodic_icdf_9_7_1d_level1(data):
+    """
+    Level 1 one-dimensional backword wavelet transformation.
+    The 0th axis is used in the transformation.
+
+    @param data (numpy.array)
+        N-dimensional array (N >= 1). This array will be *destroyed*.
+
+    @return (numpy.array)
+        Result of transformation.
+    """
+    size = len(data)
+    if size <= 1:
+        return data
+
+    assert(size % 2 == 0)
+
+    # interleave
+    scaling_size = (size + 1) // 2
+    ret = data
+    data = numpy.empty_like(ret)
+    data[0::2] = K     * ret[:scaling_size]
+    data[1::2] = (1/K) * ret[scaling_size:]
+
+    # From equations (F-7) in ITU-T Rec. T.800
+
+    # update (2)
+    data[0] -= delta * (data[1] + data[-1])
+    data[2:-1:2] -= delta * (data[1:-2:2] + data[3::2])
+
+    # predict (2)
+    data[1:-1:2] -= gamma * (data[0:-2:2] + data[2::2])
+    data[-1] -= gamma * (data[0] + data[-2])
+
+    # update (1)
+    data[0] -= beta * (data[1] + data[-1])
+    data[2:-1:2] -= beta * (data[1:-2:2] + data[3::2])
+
+    # predict (1)
+    data[1:-1:2] -= alpha * (data[0:-2:2] + data[2::2])
+    data[-1] -= alpha * (data[0] + data[-2])
+
+    return data
+
+
+def scaled_size(size, level):
+    """
+    Get the size of the level-n average portion
+    (as opposed to the detail portion).
+    @param size (int or tuple of int)
+        Size of the original vector
+    @param level (int or tuple of int)
+        Level of wavelet transformation
+    """
+    size  = numpy.asarray(size , dtype=int)
+    level = numpy.asarray(level, dtype=int)
+    return (size + (2**level - 1)) // (2**level)
diff --git a/python/lsst/obs/subaru/straylightRemoval/waveletCompression/ybackground.py b/python/lsst/obs/subaru/straylightRemoval/waveletCompression/ybackground.py
new file mode 100644
index 00000000..b8c1facb
--- /dev/null
+++ b/python/lsst/obs/subaru/straylightRemoval/waveletCompression/ybackground.py
@@ -0,0 +1,157 @@
+# Copyright (C) 2017  HSC Software Team
+# Copyright (C) 2017  Sogo Mineo
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+from __future__ import division, absolute_import, print_function
+
+from . import config
+from . import dwt
+
+import numpy
+import astropy.io.fits as pyfits
+import scipy.interpolate
+
+import os
+
+
+def main():
+    import argparse
+    parser = argparse.ArgumentParser(
+        description="""
+            Generate background image (per 1 second exposure) for a specific
+            instrument rotation angle (range).
+            We hypothesize that the instrument rotator rotates at a constant
+            angular velocity during exposure. This hypothesis has yet to be
+            verified.
+        """,
+    )
+
+    parser.add_argument("ccd", type=int, help="""
+        DET-ID for the CCD.
+    """)
+    parser.add_argument("angle_start", type=float, help="""
+        Instrument rotator angle at the start of the exposure.
+    """)
+    parser.add_argument("angle_end", nargs="?", type=float, help="""
+        Instrument rotator angle at the end of the exposure.
+        If the instrument rotator doesn't rotate during the exposure,
+        leave this argument omitted.
+    """)
+    parser.add_argument("--output", default="output.fits", help="""
+        Output file name.
+    """)
+
+    args = parser.parse_args()
+
+    img = get_ybackground(args.ccd, args.angle_start, args.angle_end)
+    pyfits.PrimaryHDU(img).writeto(args.output, clobber=True)
+
+
+def get_ybackground(ccd, angle_start, angle_end=None):
+    """
+    Get y-band background image for the CCD `ccd`.
+    It is hypothesized that the instrument rotator rotates
+    at a constant angular velocity.
+
+    @param ccd (int)
+        DET-ID of the CCD.
+    @param angle_start (float)
+        Instrument rotation angle in degrees at the start of exposure.
+    @param angle_end (float)
+        Instrument rotation angle in degrees at the end of exposure.
+
+    @return (numpy.array[][])
+    """
+    path = os.path.join(config.data_dir, "ybackground-{:03d}.fits".format(ccd))
+    hdulist = pyfits.open(path)
+
+    header = hdulist[0].header
+
+    # full-size ccd height & channel width
+    ccd_h, ch_w = header["F_NAXIS2"], header["F_NAXIS1"]
+    # saved data is compressed to 1/2**scale_level of the original size
+    image_scale_level = header["WTLEVEL2"], header["WTLEVEL1"]
+    angle_scale_level = header["WTLEVEL3"]
+
+    ccd_w = ch_w * len(hdulist)
+    ccd_img = numpy.empty(shape=(ccd_h, ccd_w), dtype=numpy.float32)
+
+    for ch, hdu in enumerate(hdulist):
+        volume = _upscale_volume(hdu.data, angle_scale_level)
+
+        if angle_end is None:
+            img = volume(angle_start)
+        else:
+            img = volume.integrate(angle_start, angle_end) * (1.0 / (angle_end - angle_start))
+
+        ccd_img[:, ch_w*ch:ch_w*(ch+1)] = _upscale_image(img, (ccd_h, ch_w), image_scale_level)
+
+    return ccd_img
+
+
+def _upscale_image(img, target_shape, level):
+    """
+    Upscale the given image to `target_shape` .
+
+    @param img (numpy.array[][])
+        Compressed image. `img.shape` must agree
+        with dwt.scaled_size(target_shape, level)
+    @param target_shape ((int, int))
+        The shape of upscaled image, which is to be returned.
+    @param level (int or tuple of int)
+        Level of multiresolution analysis (or synthesis)
+
+    @return (numpy.array[][])
+    """
+    h, w = dwt.scaled_size(target_shape, level)
+
+    large_img = numpy.zeros(shape=target_shape, dtype=float)
+    large_img[:h, :w] = img
+
+    return dwt.icdf_9_7(large_img, level)
+
+
+def _upscale_volume(volume, level):
+    """
+    Upscale the given volume (= sequence of images) along the 0-th axis,
+    and return an instance of a interpolation object that interpolates
+    the 0-th axis. The 0-th axis is the instrument rotation.
+
+    @param volume (numpy.array[][][])
+        Sequence of images.
+    @param level (int)
+        Level of multiresolution analysis along the 0-th axis.
+
+    @return (scipy.interpolate.CubicSpline)
+        You get a slice of the volume at a specific angle (in degrees)
+        by calling the returned value as `ret_value(angle)` .
+    """
+    angles = 720
+    _, h, w = volume.shape
+
+    large_volume = numpy.zeros(shape=(angles+1, h, w), dtype=float)
+
+    layers = dwt.scaled_size(angles, level)
+    large_volume[:layers] = volume
+
+    large_volume[:-1] = dwt.periodic_icdf_9_7_1d(large_volume[:-1], level, axis=0)
+    large_volume[-1] = large_volume[0]
+
+    x = numpy.arange(angles+1) / 2.0
+    return scipy.interpolate.CubicSpline(x, large_volume, axis=0, bc_type="periodic")
+
+
+if __name__ == "__main__":
+    main()
