# Order-disorder effect of $A$-site and oxygen-vacancy on magnetic and
transport properties of Y1/4Sr3/4CoO3-δ

Shun Fukushima  Tomonori Sato  Daisuke Akahoshi  Hideki Kuwahara
<EMAIL_ADDRESS>Department of Physics, Sophia University
Chiyoda-ku, Tokyo 102-8554, JAPAN

###### Abstract

We have synthesized $A$-site ordered ($A$O)- and $A$-site disordered
($A$D)-Y1/4Sr3/4CoO3-δ (YSCO) with various oxygen deficiency $\delta$, and
have made a comparative study of the structural and physical properties. We
have found that $A$-site (Y/Sr) ordering produces the unconventional oxygen-
vacancy ordered (OO) structure, and that the magnetic and transport properties
of both $A$O- and $A$D-YSCO strongly depend on the oxygen-vacancy (or excess
oxygen) ordering pattern. $A$O-YSCO with a stoichiometric $\delta$ of 0.5 has
the unconventional OO structure reflecting Y/Sr ordering pattern. With
decreasing $\delta$ from 0.5, the overall averaged OO structure is essentially
unchanged except for an increase of occupancy ratio for the oxygen-vacant
sites. At $\delta=0.34$, excess oxygen atoms are ordered to form a novel
superstructure, which is significant for the room-temperature ferromagnetism
of $A$O-YSCO. In $A$D-YSCO, on the other hand, the quite different OO
structure, which is of a brownmillerite-type, is found only in the vicinity of
$\delta=0.5$.

###### pacs:

75.30.-m, 75.30.Kz, 75.50.Dd

##  I Introduction

Transition-metal oxides with perovskite structure, $R_{1-x}Ae_{x}B$O3 ($R$:
rare-earth, $Ae$: alkaline-earth, and $B$: transition-metal), and their
derivatives exhibit rich physical properties such as high-$T_{C}$
superconductivity in Cu oxides and colossal magnetoresistance (CMR) effect in
Mn oxides.Imada_RMP_70  Perovskite-related oxides have been intensively
studied from the viewpoints of not only strongly correlated electron physics
but also potential application for correlated electron devices. In normal
perovskite structures, $R$ and $Ae$ atoms randomly occupy the $A$-sites.
Resultant $A$-site randomness often suppresses electronic phases, resulting in
a large reduction of phase transition temperatures, which makes its practical
application difficult. Therefore, it is required to reduce or remove such
$A$-site randomness for the achievement of its practical application.

Perovskite oxides with $A$-site ordered ($A$O) structures are expected as
promising candidates for correlated electron devices because they are free
from $A$-site randomness. One of the typical examples is a high-$T_{C}$
superconductor YBa2Cu3O7, which has a relatively high superconducting
transition temperature of 90 K among high-$T_{C}$ superconductors. Another
example is $A$O-$R$BaMn2O6, in which the charge- and orbital-ordering
temperature ($T_{\rm CO}=300$-480 K) is much higher than that of conventional
manganites with $A$-site disordered ($A$D) structures.Nakajima_JPCS_63 ;
Nakajima_JPSJ_71  $A$-site cation ordering not only raises the $T_{\rm CO}$
but also gives birth to an anomalous electronic phase that is not found in
$A$D perovskite manganites. In $A$O-$R$BaMn2O6, for example, the charge- and
orbital-ordering pattern and the associated magnetic structure are quite
different from those of $A$D manganites.Arima_PRB_66 ; Williams_PRB_66 ;
Uchida_JPSJ_71 ; Kageyama_JPSJ_72  Furthermore, in $A$O-$R$BaMn2O6, the
electronic phases such as the ferromagnetic metallic, charge- and orbital-
ordered, and $A$-type antiferromagnetic phases compete with each other to form
a multicritical point at $R$ = Nd. $R$/Ba disordering largely modifies the
electronic phase diagram.Akahoshi_PRL_90 ; Nakajima_JPSJ_73  In
$A$D-$R$BaMn2O6, the ferromagnetic metallic phase and charge- and orbital-
ordered one are largely suppressed, and consequently, large phase fluctuation
is enhanced near the original multicritical region ($R$ = Nd). Such large
phase fluctuation is significant for the CMR effect.Akahoshi_PRL_90 ;
Motome_PRL_91

A cobaltite with a new type of an $A$O perovskite structure, Y1-xSrxCoO3-δ has
been recently reported by Istomin et al.Istomin_Chem.Mater._15  and Withers et
al.Withers_JSSC_174  In Fig. 1(d), the crystal structure of $A$O-Y1-xSrxCoO3-δ
($x$ = 3/4, $\delta=0.50$) is displayed. Y and Sr atoms are ordered within the
$ab$-plane, and a CoO6 octahedral layer and a CoO4 tetrahedral layer
alternately stack along the $c$-axis to form four-times periodicity. Oxygen-
vacancies in the CoO4 tetrahedral layers (black spheres in Figs. 1(d) and (e))
are regularly arranged in an unconventional way: four oxygen-vacancies form a
cluster near Y atoms. The oxygen deficiency $\delta$ ($\leq 0.5$) depends on
occupancy ratio only for the oxygen-vacant sites. When the oxygen-vacant sites
are fully occupied, $\delta$ becomes zero. Kobayashi et al. reported that
$A$O-Y1-xSrxCoO3-δ with $0.75\leq x\leq 0.8$ shows a ferromagnetic behavior
above room-temperature.Kobayashi_PRB_72  Ishiwata et al. proposed that
ordering of Co $e_{g}$ orbitals causes the room-temperature
ferromagnetism.Ishiwata_PRB_75  As another characteristic of
$A$O-Y1-xSrxCoO3-δ, the physical properties are susceptible to the variation
of $\delta$. The ground state of $A$O-Y1/3Sr2/3CoO3-δ changes from an
antiferromagnetic insulator ($\delta=0.34$) to a ferromagnetic metal
($\delta=0.30$) with a slight decrease of $\delta$.Maignan_JSSC_178

In this study, we have prepared $A$O- and $A$D-Y1/4Sr3/4CoO3-δ with different
$\delta$, and have made a comparative study of their structural and physical
properties in order to reveal the effect of $A$-site and oxygen-vacancy
arrangement on the physical properties of the perovskite cobaltites.

##  II Experiment

Table 1: Annealing conditions, oxygen deficiency $\delta$ determined by
iodometric titration, and corresponding Co valences of Y1/4Sr3/4CoO3-δ (YSCO).
(a) $A$O/OO-YSCO
---
annealing condition | $\delta$ | Co valence
400 atm O2 873 K (10 h) | $0.302(17)$ | 3.15(3)
1 atm O2 773 K (2 h) | $0.340(9)$ | 3.07(2)
Ar 1173 K (12 h) | $0.437(3)$ | 2.88(1)
(b) $A$D/OD-YSCO
annealing condition | $\delta$ | Co valence
400 atm O2 873 K (10 h) | $0.152(6)$ | 3.45(1)
1 atm O2 773 K (2 h) | $0.158(4)$ | 3.43(1)
4% H2 in Ar 423 K (10 min) | $0.206(7)$ | 3.34(1)
4% H2 in Ar 473 K (10 min) | $0.222(7)$ | 3.31(1)
4% H2 in Ar 523 K (10 min) | $0.272(3)$ | 3.21(1)
4% H2 in Ar 573 K (10 min) | $0.319(5)$ | 3.11(1)
(c) $A$D/OO-YSCO
annealing condition | $\delta$ | Co valence
4% H2 in Ar 573 K (6 h) | $0.470(4)$ | 2.81(1)

$A$O- and $A$D-Y1/4Sr3/4CoO3-δ (YSCO) were prepared in a polycrystalline form
by solid state reaction. Mixed powders of Y2O3, SrCO3, and CoO were heated at
1073 K in air with a few intermediate grindings, and sintered at 1423 K in
air. The sintered powder was then treated at 1173 K in Ar to enhance $A$-site
ordering. The resultant product has the $A$O-structure with $\delta=0.44$.
$A$D-YSCO was obtained by rapidly cooling the sintered powder from 1473 (air)
to 77 K (liquid nitrogen). Hereafter, we refer to the thus obtained sample as
RC-YSCO. We controlled $\delta$ of $A$O- and $A$D-YSCO by annealing the
samples under various conditions (Table I). The values of $\delta$ were
determined by iodometric titration. Each powdered sample (ca. 30 mg) was
dissolved in 1 M HCl solution (ca. 50 ml) containing excess aqueous KI. The
amount of the formed I2 was titrated with 0.01 M Na2S2O3 solution using starch
as a colorimetric indicator. The titration was performed more than three times
for each sample. Co valences were determined by assuming that the samples have
a nominal cation ratio. The average values of $\delta$ and the corresponding
Co valences are listed in Table I. Powder X-ray diffraction (XRD) measurements
were carried out at room-temperature using a RIGAKU Rint-2100 diffractometer
with Cu $K\alpha$ radiation. The crystal structures of the samples were
analyzed by Rietveld method using RIETAN-2000.rietan  Magnetic and transport
properties were measured from 5 to 350 K using a Quantum Design, Physical
Property Measurement System (PPMS). Magnetization measurements above room-
temperature ($T=300$-800 K) were carried out using a Quantum Design, Magnetic
Property Measurement System (MPMS).

##  III Results

###  III.1 Crystal structures

Figure 1: (Color online) X-ray-diffraction patterns of (a) $A$-site and
oxygen-vacancy ordered ($A$O/OO)-Y1/4Sr3/4CoO2.66, (b) $A$-site and oxygen-
vacancy disordered ($A$D/OD)-Y1/4Sr3/4CoO2.84, and (c) $A$-site disordered and
oxygen-vacancy ordered ($A$D/OO)-Y1/4Sr3/4CoO2.53. Solid circles and solid
lines represent observed and calculated diffraction profiles, respectively.
Vertical marks under the diffraction profiles indicate calculated peak
positions. The insets show the diffraction profiles at low angles. For
comparison, diffraction profiles of $\delta=0.30$ and 0.44 ($A$O/OO) are also
depicted in the inset of (a). Crystal structure and its idealized oxygen-
vacant layer of (d),(e) $A$O/OO-Y1/4Sr3/4CoO2.5 and (f),(g)
$A$D/OO-Y1/4Sr3/4CoO2.5. Small black spheres of (d)-(g) represent oxygen-
vacancies.

Figure 1(a) shows the XRD pattern of $A$O-YSCO with $\delta=0.34$. It should
be noted that the Bragg peak around $2\theta=11.5$ deg. is a piece of evidence
that both Y/Sr and oxygen-vacancies are ordered. Except for an additional
Bragg peak indicated by an arrow in the inset of Fig. 1(a), the XRD patterns
of $A$O-YSCO with $\delta=0.30$ and 0.44 are similar to that of $\delta=0.34$
($A$O/OO), and all the XRD patterns can be fitted to the $A$O/OO structure
($I4/mmm$) shown in Fig. 1(d). Note that all $A$O/OO-YSCO prepared in this
study, whose $\delta$ is smaller than the stoichiometric composition
($\delta=0.5$), have excess oxygen atoms. The result of the XRD measurement
indicates that these excess oxygen atoms randomly occupy the oxygen-vacant
sites (or raise the occupancy ratio of the site) in the CoO4 layers with
keeping the overall averaged OO structure shown in Figs. 1(d) and (e), that
is, the OO structure is robust against the variation of $\delta$. The
additional Bragg peak observed in $\delta=0.34$ ($A$O/OO) can be indexed as
(1/4 1/4 0) in an $a_{p}\times a_{p}\times a_{p}$ setting ($a_{p}$ denotes a
pseudo-cubic perovskite cell), evidencing the existence of a four-times
superstructure along the [1 1 0] direction. Because the additional peak is
sensitive to $\delta$ and found only in the XRD pattern of $\delta=0.34$
($A$O/OO), it is reasonable to conclude that the superstructure arises from
ordering of the excess oxygen atoms occupying the oxygen-vacant sites. This
means that $A$O-YSCO probably has another oxygen stoichiometry near
$\delta=0.34$ besides the stoichiometric $\delta=0.5$. $\delta=1/3$ (or
$3/8$), which is close to 0.34, is likely to be the second stoichiometric
composition. A similar superlattice peak (1/4 1/4 0) is also observed in
$A$O-Er0.78Sr0.22CoO2.63 by use of a high-intensity synchrotron X-ray source
Ishiwata_PRB_75 , implying that $\delta=0.34$ ($A$O/OO) has the similar oxygen
superstructure to $A$O-Er0.78Sr0.22CoO2.63. The intensity of the superlattice
peak of our sample is much stronger than that of
$A$O-Er0.78Sr0.22CoO2.63fnote_ESCO , indicating that the oxygen deficiency
$\delta$ of our sample is closer to the stoichiometry than that of
$A$O-Er0.78Sr0.22CoO2.63. The detailed crystal structure of $\delta=0.34$
($A$O/OO) is now under investigation.

Then, we exhibit the XRD pattern of RC-YSCO with $\delta=0.16$ in Fig. 1(b).
What is a significant difference from the XRD patterns of $A$O-YSCO system is
that the characteristic Bragg peaks arising from $A$-site and oxygen-vacancy
ordering totally vanish. The XRD pattern can be well fitted to a simple cubic
perovskite structure ($Pm\bar{3}m$), indicating that Y and Sr atoms randomly
occupy the $A$-sites, and that oxygen-vacancies are randomly distributed in
the structure; RC-YSCO has an $A$D and oxygen-vacancy disordered (OD)
structure. All RC-YSCO with $0.15\leq\delta\leq 0.32$ have the same simple
cubic structure as RC-YSCO with $\delta=0.16$. From now on, we will refer to
these compounds as $A$D/OD-YSCO.

On the other hand, the XRD pattern of RC-YSCO with $\delta$ = 0.47 (Fig. 1(c))
is quite different from those of RC-YSCO with $0.15\leq\delta\leq 0.32$. A
characteristic Bragg peak, which implies $A$-site and/or oxygen-vacancy
ordering, is clearly seen around $2\theta=11.5$ deg. By annealing RC-YSCO with
$\delta=0.47$ in O2 at 573 K, which is low enough to prevent $A$-site cations
from moving around to rearrange, the characteristic peak totally disappears,
and the annealed compound takes a simple cubic perovskite structure.
Therefore, we conclude that only oxygen-vacancy ordering contributes to the
appearance of the characteristic peak; RC-YSCO with $\delta=0.47$ has an
$A$D/OO structure. Rietveld analysis reveals that the OO structure of
$A$D-YSCO with $\delta=0.47$ is of a brownmillerite (Ca2(Fe,Al)2O5)-type (Fig.
1(f) and (g), $Ibm2$)Colville_Acta_Cryst_B27 , which is quite different from
that of $A$O/OO-YSCO shown in Figs. 1(d) and (e).

To summarize this section, YSCO can be classified into three types of the
crystal structures: the $A$O/OO, $A$D/OD, and $A$D/OO structures. In the
following section, we demonstrate that the arrangement of Y/Sr and oxygen-
vacancies considerably affects the physical properties of YSCO.

###  III.2 $A$-site and oxygen-vacancy ordered ($A$O/OO)-YSCO

Figure 2: (Color online) Temperature dependence of (a) magnetization, (b)
resistivity, and (c) inverse magnetization ($H/M$) of $A$O/OO-Y1/4Sr3/4CoO3-δ.
FC: field cooled.

Figures 2(a) and (b) show the temperature dependence of the magnetization and
resistivity of $A$O/OO-YSCO with $\delta=0.30$, 0.34, and 0.44. The
magnetization of $\delta$ = 0.34 ($A$O/OO) shows a weak ferromagnetic
transition around 330 K as previously reported by Kobayashi et
al.Kobayashi_PRB_72  Ishiwata et al. suggest that the weak ferromagnetism (the
room-temperature ferromagnetism) originates in canted antiferromagnetic or
ferrimagnetic order caused by ordering of Co $e_{g}$ orbital.Ishiwata_PRB_75
Then, the magnetization of $\delta=0.34$ ($A$O/OO) shows a sharp cusp around
300 K, implying that a magnetic or spin-state transition occurs. The
resistivity of $\delta=0.34$ ($A$O/OO) exhibits an anomaly around the weak
ferromagnetic transition temperature, below which it is insulating. Figure
2(c) shows the inverse magnetization of $A$O/OO-YSCO with $\delta=0.30$, 0.34,
and 0.44. All the inverse magnetizations approximately obey the Curie-Weiss
law above 400 K. The Weiss temperature $\theta_{W}$ and effective moment
$P_{\rm eff}$ of $\delta=0.34$ ($A$O/OO) are found to be $-60$ K and 3.3
$\mu_{\rm B}/{\rm Co}$, respectively. The magnetic and transport properties of
$\delta=0.34$ ($A$O/OO) observed in this study are consistent with those
previously reported by Kobayashi et al.Kobayashi_PRB_72 ; Kobayashi_JPSJ_75 ,
except for the magnetic cusp at 300 K. The origin of the cusp will be
discussed later.

Figure 3: (Color online) Magnetic field dependence of magnetization ($M$-$H$
curves) of $A$O/OO-Y1/4Sr3/4CoO3-δ at 5 K.

In $\delta=0.30$ ($A$O/OO), the weak ferromagnetic insulating phase above
room-temperature is considerably suppressed, while the magnetization is
continuously increasing with decreasing temperature, and the ferromagnetic
correlation is larger than that of $\delta=0.34$ ($A$O/OO) below $\sim 100$ K
(Figs. 2(a) and 3). With decreasing $\delta$ from 0.34 to 0.30, the
$\theta_{W}$ increases from $-60$ K to 60 K, that is, an antiferromagnetic
interaction between Co ions turns into a ferromagnetic one. The resistivity of
$\delta=0.30$ ($A$O/OO) largely drops in the whole temperature region compared
with that of $\delta=0.34$ ($A$O/OO). These results show that $A$O-YSCO
approaches a ferromagnetic metal with a decrease of $\delta$, i.e. with an
increase of Co valence, and that the weak ferromagnetic insulating phase
($\delta\approx 0.34$) and ferromagnetic metallic clusters (oxygen rich
region: $\delta\ll 0.34$) coexist in $\delta=0.30$ ($A$O/OO).

On the other hand, the magnetization of $\delta=0.44$ ($A$O/OO), which is
close to the stoichiometric composition of $\delta=0.5$ and has excess oxygen
atoms occupying 12 % of the oxygen-vacant sites, shows a slight increase
around 230 K, suggesting that a weak ferromagnetic transition occurs. However,
as seen from the magnetic field dependence of the magnetization ($M$-$H$
curves) (Fig. 3), the weak ferromagnetic magnetization is much smaller than
that of $\delta=0.34$ ($A$O/OO). The inverse magnetization of $\delta=0.44$
($A$O/OO) shows a clear anomaly around 350 K (Fig. 2(c)), which is attributed
to the remnant of the room-temperature ferromagnetic phase most stabilized
around $\delta=0.34$. These results indicate that the magnetic properties of
$\delta=0.44$ ($A$O/OO) can be explained by the coexistence of the matrix
phase ($\delta=0.5$) and embedded clusters ($\delta\approx 0.34$) due to the
excess oxygen atoms. The weak ferromagnetism below 230 K probably comes from
the matrix phase.

To summarize this section, the room-temperature ferromagnetism is observed
only in $\delta=0.34$ ($A$O/OO), and a slight deviation of $\delta$ from 0.34
strongly suppresses the room-temperature ferromagnetism. The sensitivity of
the room-temperature ferromagnetic phase to $\delta$ also supports our
aforementioned conclusion that $A$O/OO-YSCO has another oxygen stoichiometry
around $\delta=0.34$.

###  III.3 $A$-site and oxygen-vacancy disordered ($A$D/OD)-YSCO

Figure 4: (Color online) Temperature dependence of (a) magnetization and (b)
resistivity of $A$D/OD-Y1/4Sr3/4CoO3-δ. The inset of (a) shows the inverse
magnetization of $A$D/OD-Y1/4Sr3/4CoO3-δ ($\delta=0.15$)

We display in Figs. 4(a) and (b) the temperature dependence of the
magnetization and resistivity of $A$D/OD-YSCO with $0.15\leq\delta\leq 0.32$.
In $\delta=0.15$ ($A$D/OD), the resistivity exhibits metallic transport,
except for a slight upturn at low temperatures, and the magnetization abruptly
increases below 160 K. The $\theta_{W}$ ($=190$ K) estimated from the inverse
magnetization (the inset of Fig. 4(a)) is close to the magnetic transition
temperature. The $M$-$H$ curve at 5 K exhibits a ferromagnetic behavior with a
saturation magnetization of 1.2 $\mu_{\rm B}/{\rm Co}$ (Fig. 5). These facts
indicate that $\delta=0.15$ ($A$D/OD) is a typical ferromagnetic metal. Such a
ferromagnetic metallic behavior is sometimes observed in conventional
perovskite cobaltites with $A$D structure such as La1-xSrxCoO3 ($x\geq
0.18$).Itoh_JPSJ_63 ; Senaris_JSSC_118  It is widely accepted that the origin
of the ferromagnetic metallicity of La1-xSrxCoO3 can be explained by the
double exchange interaction.Itoh_JPSJ_63 ; Kriener_PRB_69 ; Saitoh_PRB_56 ;
Yamaguchi_JPSJ_64  The saturation magnetization of $\delta=0.15$ ($A$D/OD) is
close to that of La0.8Sr0.2CoO3 ($\sim 1.3$ $\mu_{\rm B}/{\rm Co}$ at 5 K),
indicating that Co3+ and Co4+ of $\delta=0.15$ ($A$D/OD) in the ferromagnetic
metallic state have similar electronic configurations to those of
La1-xSrxCoO3.Senaris_JSSC_118 ; Kriener_PRB_69  Thus, it is reasonable to
conclude that the ferromagnetic metallic state of $A$D-YSCO with $\delta=0.15$
is stabilized via the double-exchange mechanism.

Figure 5: (Color online) Magnetic field dependence of magnetization ($M$-$H$
curves) of $A$D/OD-Y1/4Sr3/4CoO3-δ at 5 K.

With an increase of $\delta$, the ferromagnetic metallic state is steeply
suppressed (Fig. 4). In $\delta=0.27$ and 0.32 ($A$D/OD), the ferromagnetic
component is negligible as seen from the temperature dependence of the
magnetization and the $M$-$H$ curves, and the resistivity exhibits a typical
insulating behavior (Figs. 4 and 5).

###  III.4 $A$-site disordered and oxygen-vacancy ordered ($A$D/OO)-YSCO

Figure 6(a) shows the temperature dependence of the magnetization and
resistivity of $A$D/OO-YSCO with $\delta=0.47$, which has the brownmillerite-
type OO structure as described in section III.1. The magnetization of $A$D/OO-
YSCO with $\delta=0.47$ abruptly increases around 130 K just like that of
$A$D/OD-YSCO with $\delta=0.15$ (Fig. 4(a)). However, the saturation
magnetization of $0.15$ $\mu_{\rm B}/{\rm Co}$ (Fig. 6(b)) is much smaller
than that of $A$D/OD-YSCO with $\delta=0.15$, and the resistivity of $A$D/OO-
YSCO with $\delta=0.47$ shows an insulating behavior in the whole temperature
region. The parent compound of $A$D/OO-YSCO, Sr2Co2O5, which also has the
brownmillerite-type OO structure, undergoes a similar weak ferromagnetic
transition at 200 K and a $G$-type antiferromagnetic transition at 537
K.Munoz_PRB_78  In $A$D/OO-YSCO with $\delta=0.47$, similarly, a $G$-type
anitiferromagnetic transition might occur far above room-
temperature.fnote_HiTemp  We note that the weak ferromagnetic moment of
$A$D/OO-YSCO with $\delta=0.47$ is more than 10 times larger than that of
Sr2Co2O5. In $A$D/OO-YSCO with $\delta=0.47$, which is in a mixed valence
state of Co3+/Co2+, a ferrimagnetic transition accompanying charge-order might
occur at 130 K.

Figure 6: (a) Temperature dependence of magnetization and resistivity, and (b)
magnetic field dependence of magnetization ($M$-$H$ curve) of
$A$D/OO-Y1/4Sr3/4CoO3-δ ($\delta=0.47$).

##  IV Discussion

Now we discuss the effect of the $A$-site and oxygen-vacancy arrangement on
the physical properties of YSCO. $A$O-YSCO with the stoichiometric composition
of $\delta=0.5$ has the unconventional OO structure, in which, four oxygen-
vacancies in the oxygen deficient CoO2-2δ layers form a cluster near Y atom in
the adjacent Y1/4Sr3/4O layers as shown in Figs. 1(d) and (e). Y/Sr
disordering largely modifies the arrangement of the oxygen-vacancies.
$A$D-YSCO with $\delta=0.5$ has the brownmillerite-type OO structure (Figs.
1(f) and (g)), which is often found in conventional perovskite oxides with
large oxygen deficiency such as Sr2Fe2O5Harder_BM ,
Ca2Fe2O5Colville_Acta_Cryst_B26 , and Sr2Co2O5Takeda_JPSJ_33 . From these
results, it is obvious that Y/Sr ordering gives rise to the unconventional OO
structure reflecting the Y/Sr ordering pattern. The OO structure of $A$O-YSCO
is so robust that excess oxygen atoms partially occupy the oxygen-vacant sites
with keeping the overall OO structure shown in Figs. 1(d) and (e). In
contrast, the OO structure of $A$D-YSCO (Figs. 1(f) and (g)) is so fragile
that it is easily destroyed by introduction of a small amount of excess oxygen
atoms. Consequently, the OO structure is found only in the vicinity of
$\delta=0.5$ in the case of $A$D-YSCO. The robustness of the OO structure of
$A$O-YSCO originates in periodic potential due to $A$-site ordering.

As for the magnetic properties, both $A$O/OO- and $A$D/OO-YSCO with the
stoichiometric $\delta$ (around 0.34 and 0.47 respectively) exhibit the weak
ferromagnetic behaviors as demonstrated in section III. It should be noted
that weak ferromagnetic behaviors are also observed in other oxygen
stoichiometric $A$O perovskites such as YBaCo2O5+x ($x=0.50$ and 0.44) and
YBaMn2O5+x ($x=0.50$), which have the superstructures formed by excess oxygen
atoms.Akahoshi_JSSC_156 ; Karppinen_JSSC_177  Therefore, it is probable that
oxygen-vacancy (or excess oxygen) ordering plays a significant role in the
weak ferromagnetism of both $A$O/OO- and $A$D/OO-YSCO.

Then, we propose one plausible model to explain the origin of the weak
ferromagnetism. In the following discussion, we suppose that the magnetic
interaction between Co ions is antiferromagnetic. Considering a CoO6
octahedron adjacent to a CoO4 tetrahedron (or a CoO5 square pyramid) in oxygen
deficient perovskite structures, an inversion symmetry is absent at the center
between the two Co ions. In such a case, a spin-canted moment is locally
induced through the Dzyaloshinskii-Moriya interaction.Dzyaloshinskii_JPCS_4 ;
Moriya_PhysRev_120  Another possible case is as follows: if the spin states of
CoO6 and CoO4 (or CoO5) were different from each other, the CoO6-CoO4 (or
CoO5) pair could have a local ferrimagnetic moment. In the OO structures, the
CoO6-CoO4 (or CoO5) pairs, i.e., the local net moments are regularly arranged.
As a result, the weak ferromagnetic magnetization is macroscopically observed.
On the other hand, in the OO structure with large oxygen nonstoichiometry
(excess oxygen randomly occupy the oxygen-vacant sites) or in the OD
structures ($\delta=0.27$ and 0.32 ($A$D/OD)) in which, the CoO6-CoO4 (or
CoO5) pairs are randomly distributed, the local net moments oriented randomly
are canceled out with each other. The reason why the weak ferromagnetic moment
of $\delta=0.44$ ($A$O/OO) below 230 K is very small may be attributed to
randomness due to excess oxygen atoms occupying 12 % of the oxygen-vacant
sites. Note that the ferromagnetism of $A$D/OD-YSCO with $\delta=0.15$ is
induced through the double-exchange interaction as described in section III.3.

The magnetization of $\delta=0.34$ ($A$O/OO) we prepared exhibits the sharp
cusp around 300 K (Fig. 2(a)), which is not clearly discerned in $A$O-YSCO
previously reported.Kobayashi_PRB_72  In this study, $A$O-YSCO was treated in
Ar to enhance Y/Sr ordering as mentioned in II, while $A$O-YSCO in the
previous reports was not. As a result, the degree of Y/Sr order of our sample
is higher than that of $A$O-YSCO in the previous reports, which may be the
main cause of the emergence of the cusp. Another plausible explanation for the
emergence of the cusp is that our sample is closer to the oxygen stoichiometry
than $A$O-YSCO in the previous reports.

##  V Summary

We have synthesized $A$-site ordered ($A$O)- and $A$-site disordered
($A$D)-Y1/4Sr3/4CoO3-δ (YSCO) with various oxygen deficiency $\delta$, and
have made a comparative study of their structural and physical properties.
$A$O-YSCO with $\delta=0.44$, which is near the stoichiometry of $\delta=0.5$,
has the unconventional oxygen-vacancy ordered (OO) structure reflecting
$A$-site ordering pattern, and shows the weak ferromagnetic behavior below 230
K. Decreasing $\delta$ from 0.5 inducing an increase in the occupancy ratio of
the oxygen-vacant site does not essentially change the overall OO structure.
In $A$O-YSCO with $\delta=0.34$, which is close to another oxygen
stoichiometry, excess oxygen ordering in the oxygen-vacant sites causes the
room-temperature ferromagnetism. On the other hand, in $A$D-YSCO, the
different OO structure, which is of the brownmillerite-type, is found only in
the vicinity of $\delta=0.5$, and also shows the weak ferromagnetic behavior
below 130 K. $A$D/OD-YSCO with $\delta=0.15$ exhibits a ferromagnetic metallic
behavior which is attributed to the double-exchange interaction. With an
increase of $\delta$, that is, with a decrease of Co valence, the
ferromagnetic metallic phase is suppressed. Our main conclusions are as
follows: First, $A$-site ordering gives rise to the unconventional OO
structure, and makes the OO structure robust against randomness due to excess
oxygen atoms. Second, oxygen-vacancy (or excess oxygen) ordering is
indispensable for the weak ferromagnetic behaviors of both $A$O- and
$A$D-YSCO.

###### Acknowledgements.

We thank D. Vieweg and A. Loidl for discussions and their help with the
magnetization measurements using MPMS. This work was partially supported by
Iketani Science and Technology Foundation, the Mazda Foundation, the Asahi
Glass Foundation, and by Grant-in-Aid for Scientific Research (C) from Japan
Society for the Promotion of Science (JSPS).

## References

  * (1) M. Imada, A. Fujimori, and Y. Tokura, Rev. Mod. Phys. 70, 1039 (1998).
  * (2) T. Nakajima, H. Kageyama, and Y. Ueda, J. Phys. Chem. Solids 63, 913 (2002).
  * (3) T. Nakajima, H. Kageyama, H. Yoshizawa, and Y. Ueda, J. Phys. Soc. Jpn. 71, 2843 (2002).
  * (4) T. Arima, D. Akahoshi, K. Oikawa, T. Kamiyama, M. Uchida, Y. Matsui, and Y. Tokura, Phys. Rev. B 66, 140408(R) (2002).
  * (5) A. J. Williams and J. P. Attfield, Phys. Rev. B 66, 220405(R) (2002).
  * (6) M. Uchida, D. Akahoshi, R. Kumai, Y. Tomioka, T. Arima, Y. Tokura, and Y. Matsui, J. Phys. Soc. Jpn. 71, 2605 (2002).
  * (7) H. Kageyama, T. Nakajima, M. Ichihara, Y. Ueda, H. Yoshizawa, and K. Ohoyama, J. Phys. Soc. Jpn. 72, 241 (2003).
  * (8) D. Akahoshi, M. Uchida, Y. Tomioka, T. Arima, Y. Matsui, and Y. Tokura, Phys. Rev. Lett. 90, 177203 (2003).
  * (9) T. Nakajima, H. Yoshizawa, and Y. Ueda, J. Phys. Soc. Jpn. 73, 2283 (2004).
  * (10) Y. Motome, N. Furukawa, and N. Nagaosa, Phys. Rev. Lett. 91, 167204 (2003).
  * (11) S. Ya. Istomin, J. Grins, G. Svensson, O. A. Drozhzhin, V. L. Kozhevnikov, E. V. Antipov, and J. P. Attfield, Chem. Mater. 15, 4012 (2003).
  * (12) R. L. Withers, M. James, and D. J. Goossens, J. Solid State Chem. 174, 198 (2003).
  * (13) W. Kobayashi, S. Ishiwata, I. Terasaki, M. Takano, I. Grigoraviciute, H. Yamauchi, and M. Karppinen, Phys. Rev. B 72, 104408 (2005).
  * (14) S. Ishiwata, W. Kobayashi, I. Terasaki, K. Kato, and M. Takata, Phys. Rev. B 75, 220406(R) (2007).
  * (15) A. Maignan, S. Hébert, V. Caignaert, V. Pralong, and D. Pelloquin, J. Solid State Chem. 178, 868 (2005).
  * (16) F. Izumi and T. Ikeda, Mater. Sci. Forum 198, 321 (2000).
  * (17) Ishiwata et al. stated that the superstructure peak could not have been detected without a high-intensity synchrotron X-ray source (Ref. Ishiwata_PRB_75, ), while we could detect it by use of a commercial X-ray source.
  * (18) A. A. Colville and S. Geller, Acta. Crystallogr. Sect. B 27, 2311 (1971).
  * (19) W. Kobayashi, S. Yoshida, and I. Terasaki, J. Phys. Soc. Jpn. 75, 103702 (2006).
  * (20) M. Itoh, I. Natori, S. Kubota, and K. Motoya, J. Phys. Soc. Jpn. 63, 1486 (1994).
  * (21) M. A. Señarís-Rodríguez and J. B. Goodenough, J. Solid State Chem. 118, 323 (1995).
  * (22) M. Kriener, C. Zobel, A. Reichl, J. Baier, M. Cwik, K. Berggold, H. Kierspel, O. Zabara, A. Freimuth, and T. Lorenz, Phys. Rev. B 69, 094417 (2004).
  * (23) S. Yamaguchi, H. Taniguchi, H. Takagi, T. Arima, and Y. Tokura, J. Phys. Soc. Jpn. 64, 1885 (1995).
  * (24) T. Saitoh, T. Mizokawa, A. Fujimori, M. Abbate, Y. Takeda, and M. Takano, Phys. Rev. B 56, 1290 (1997).
  * (25) A. Muñoz, C. de la Calle, J. A. Alonso, P. M. Botta, V. Pardo, D. Baldomir, and J. Rivas, Phys. Rev. B 78, 054404 (2008).
  * (26) Since $A$D/OO-YSCO is chemically unstable, we could not obtain reliable data of the physical properties at higher temperatures.
  * (27) M. Harder and H. Müller-Buschbaum, Z. Anorg. Allg. Chem. 464, 169 (1980).
  * (28) A. A. Colville, Acta. Crystallogr. Sect. B 26, 1469 (1970).
  * (29) T. Takeda, Y. Yamaguchi, and H. Watanabe, J. Phys. Soc. Jpn. 33, 970 (1972).
  * (30) D. Akahoshi and Y. Ueda, J. Solid State Chem. 156, 355 (2001).
  * (31) M. Karppinen, H. Okamoto, H. Fjellvåg, T. Motohashi, and H. Yamauchi, J. Solid State Chem. 177, 2122 (2004).
  * (32) I. Dzyaloshinskii, J. Phys. Chem. Solids 4, 241 (1958).
  * (33) T. Moriya, Phys. Rev. 120, 91 (1960).


# Image Denoising in FPGA using Generic Risk Estimation

Rinson Varghese  College of Engineering and Management, Punnapra
Kerala, India
Email<EMAIL_ADDRESS>Chandrasekhar Seelamantula  IISc,
Bangalore
India
Email<EMAIL_ADDRESS>Rathna G N  IISc, Bangalore
India
Email<EMAIL_ADDRESS>Ashutosh Gupta  SAC, ISRO
India
Email<EMAIL_ADDRESS>Debajyoti Dhar  SAC, ISRO
India
<EMAIL_ADDRESS>
###### Abstract

The generic risk estimator addresses the problem of denoising images corrupted
by additive white noise without placing any restriction on the statistical
distribution of the noise. In this paper, we discuss an efficient FPGA
implementation of this algorithm. We use the undecimated Haar wavelet
transform with shrinkage parameters for each sub-band as the denoising
function. The computational complexity and memory requirement of the algorithm
is first analyzed. To optimize the performance, a combination of convolution
and recursion is employed to realize Haar filter bank and gradient descent
algorithm is used to find the shrinkage parameters. A fully pipelined and
parallel architecture is developed to achieve high throughput. The proposed
design achieves an execution time of 3.5ms for an image of size $512\times
512$. We also show that the recursive implementation of Haar wavelet is more
expensive than the direct implementation in terms of hardware utilization.

##  I Introduction

Images and videos are omnipresent in our lives. We see them on television,
surveillance systems, medical imaging systems etc. In every image processing
application, we extract some information from one or more images taken by the
camera and process it further. Unfortunately, the fundamental physical laws
which govern the conversion of light into pixel values in Charge Coupled
Device (CCD) sensors introduce noise into the image [1]. There are several
potential noise sources in a CCD camera, such as shot noise, dark current,
pixel non-uniformity and quantization noise. In very low light applications,
when there is a time constraint on exposure time, CCD senors are seriously
degraded by their readout noise [2]. Hence the level of noise depends not only
on the sensor, but on the ambient conditions as well. Estimation and removal
of noise is thus essential to increase the performance of image processing
algorithms in the subsequent stages.

Since noise is random in nature, it is characterized in the statistical
setting described above by its distribution function. Additive White Gaussian
noise (AWGN) is the most commonly used model in image processing. In fact, the
dark current, which is dependent on the temperature and exposure time is
modeled as a Gaussian distribution. Noise may follow other distributions as
well, for instance, Shot and Read noise caused by on chip amplifiers follow
Poisson distribution where as hot pixels follow impulsive distribution [3].
Noise in magnetic resonance imaging is found to be Rician distributed [4] and
in SAR images it is speckle noise [5]. Hence there is a need for a denoising
method which is independent of the noise distribution. Further, to use these
algorithms in real time applications, we may have to meet some timing and
power constraints. Hence we go for hardware implementation of the algorithm
instead of software. FPGA is a good choice in this context since we can employ
pipelining and parallelism to improve the speed at a reduced power.

There are different works which show the application of FPGA in image
denoising and image processing. In [6], some simple filters like median
filter, smoothing filter and edge detection filter are implemented in FPGA.
Different architectures of a bilateral filter are discussed in [7, 8, 9],
where [9] is the first scalable FPGA implementation of the bilateral filter. A
reconfigurable gaussian noise detection in images is proposed in [10], where
the real time dynamic reconfiguration of the algorithm is achieved through the
processor core. In [11], an adaptive image denoising IP core based on FPGA is
proposed where adaptive gaussian filter is employed for noise elimination. The
denoising framework in most of these algorithms puts some restriction on the
noise distribution. With the popularity of the wavelet transform for the past
three decades, several algorithms have been developed in wavelet domain. Only
a few of these algorithms have been implemented in FPGA. A hardware
architecture for wavelet domain video denoising system is proposed in [12]
where the denoising is performed in wavelet domain using spatially adaptive
bayesian shrinkage.

In this work we develop a hardware architecture for generic risk estimator
(GenRE) proposed in [13], where there is no restriction assumed on the noise
distribution. We use undecimated Haar wavelet transform together with
shrinkage parameters as the denoising function.

Figure 1: A block diagram of an undecimated wavelet based denoising for two
levels of decomposition with shrinkage parameters $\alpha_{i}$.

##  II Generic Risk Estimator

The objective of any image denoising algorithm is to estimate the unknown true
image from a noisy observation. We model this problem as follows.

$\mathbf{y}=\mathbf{x}+\mathbf{w}$ (1)

where,
$\mathbf{x}\in\mathbb{R^{N}}$ is the column-vectorized noise free image,
$\mathbf{w}\in\mathbb{R^{N}}$ is the noise, whose first and second order
statistics is known and assumed to have zero mean and variance.
$\mathbf{y}\in\mathbb{R^{N}}$ is the noisy observation.
We also assume that $\mathbf{w}$ is independent of $\mathbf{x}$ and the
distribution of $\mathbf{w}$ is unknown. The goal is to find an optimal
denoising function $\mathbf{f}:\mathbb{R^{N}}\mapsto\mathbb{R^{N}}$, which
gives $\hat{\mathbf{x}}$, an estimate of the true image , when operated on
$\mathbf{y}$. The optimality of the denoising function is determined by mean
square error (MSE). Hence the whole problem can be restated as follows: Find a
denoising function $\mathbf{f}$ such that the MSE

$\mathbf{C}(\mathbf{f})=\frac{1}{N}\mathbb{E}(||\mathbf{x}-\hat{\mathbf{x}}||^{2})$
(2)

between $\mathbf{x}$ and $\hat{\mathbf{x}}=\mathbf{f}(\mathbf{y})$ is
minimized. However, non-availability of ground truth becomes an obstacle for
direct calculation of the MSE. Hence instead of minimizing MSE, we find an
unbiased estimate of the MSE and minimize it to find the denoising function.

######  Theorem 1.

Let $\mathbf{f}$ be a linear denoising function, given as
$\mathbf{f}(\mathbf{y})=\mathbf{B}\mathbf{y}$,
$\mathbf{B}\in\mathbb{R}^{N{\times}N}$, then

$\text{$\mathbf{GenRE}$}(\mathbf{f})=\frac{1}{N}(||\mathbf{f}(\mathbf{y})-\mathbf{y}||^{2}+2\sigma^{2}\text{Trace}(\mathbf{B})-N\sigma^{2})$
(3)

is an unbiased estimator of the MSE.

A formal proof of the theorem is given in [13].

Figure 2: A block diagram of GenRE based denoising using UWT.

##  III GenRE based denoising using undecimated wavelet transform

In discrete wavelet transform (DWT), at every level of decomposition we down-
sample the signal which removes redundant information. Similar to this we up-
sample the signal during reconstruction which restores the original sampling
rate. Undecimated wavelet transform (UWT) removes this sampling rate
conversion and instead of down-sampling the signal, the filter itself is up-
sampled. The filter coefficients are upsampled by a factor of $2^{j-1}$ in the
$j^{\text{th}}$ level of the algorithm. Fig. 1 illustrates an undecimated Haar
wavelet based denoising which uses two levels of decomposition. The denoising
function is a combination of analysis filter bank, transform-domain shrinkage,
and synthesis filter bank. In matrix form we represent the analysis filter
bank as
$\mathbf{D}=[\mathbf{D}_{1}^{\text{T}},\mathbf{D}_{2}^{\text{T}},\dots\mathbf{D}_{j}^{\text{T}}]^{\text{T}}$
and synthesis filter bank as
$\mathbf{R}=[\mathbf{R}_{\text{1}},\mathbf{R}_{\text{2}},\dots\mathbf{R}_{\text{j}}]$,
such that $\mathbf{R}\mathbf{D}=\mathbf{I}$. Here both $\mathbf{R}_{\text{i}}$
and $\mathbf{D}_{\text{i}}$ are circulant matrices. We can describe the
overall transform domain process as follows

$\displaystyle\hat{\mathbf{x}}$ $\displaystyle=$
$\displaystyle\mathbf{f}(\mathbf{y})=\sum_{i=1}^{j}\mathbf{R}_{i}\boldsymbol{\alpha}_{i}\mathbf{D}_{i}\mathbf{y}=\sum_{i=1}^{j}\boldsymbol{\alpha}_{i}\mathbf{H}_{i}\mathbf{y}=\mathbf{B}\mathbf{y}$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{j}\boldsymbol{\alpha}_{i}\mathbf{\psi}_{i}=\mathbf{\Psi}\boldsymbol{\alpha}$
(4)

where $\mathbf{H}_{i}=\mathbf{R}_{i}\mathbf{D}_{i}$, $\alpha_{i}$ is the
shrinkage parameter for $\text{i}^{\text{th}}$ sub-band and
$\\{\mathbf{\psi}_{i}\\}_{i=1}^{J+1}$ forms the columns of $\mathbf{\Psi}$.
Here

$\mathbf{B}=\sum_{\text{i =
1}}^{\text{j}}\alpha_{i}\mathbf{H}_{i}=\sum_{\text{i =
1}}^{\text{j}}\alpha_{i}\mathbf{R}_{i}\mathbf{D}_{i}$ (5)

$\mathbf{B}$ is a circulant matrix, since it is a product of two circulant
matrices. Now we can rewrite (3) as

$\displaystyle\text{$\mathbf{GenRE}$}(\mathbf{f})$ $\displaystyle=$
$\displaystyle\frac{1}{N}(||\mathbf{f}(\mathbf{y})-\mathbf{y}||^{2}+2\sigma^{2}\text{Trace}(\mathbf{B})-N\sigma^{2})$
(6) $\displaystyle=$
$\displaystyle\frac{1}{N}(||\mathbf{\Psi}\boldsymbol{\alpha}-\mathbf{y}||^{2}+2\sigma^{2}N\mathbf{B}_{11}-N\sigma^{2})$
$\displaystyle=$
$\displaystyle\frac{1}{N}(||\mathbf{\Psi}\boldsymbol{\alpha}-\mathbf{y}||^{2}+2\sum_{\text{i
= 1}}^{\text{j}}\alpha_{i}\mathbf{H}_{11}^{i}N\sigma^{2}-N\sigma^{2})$
$\displaystyle=$
$\displaystyle\frac{1}{N}(||\mathbf{\Psi}\boldsymbol{\alpha}-\mathbf{y}||^{2}+2\boldsymbol{\alpha}^{\text{T}}\mathbf{q}-N\sigma^{2})$

We can find the shrinkage parameter $\boldsymbol{\alpha}$ by minimizing
$\mathbf{GenRE}$, which gives

$\nabla_{\alpha}(\text{$\mathbf{GenRE}$}(\mathbf{f}))=2\mathbf{\Psi}^{T}\mathbf{\Psi}\boldsymbol{\alpha}^{*}-2\boldsymbol{\alpha}^{\text{T}}\mathbf{y}-2\mathbf{q}=0$
(7)

and

$\boldsymbol{\alpha}^{*}=(\mathbf{\Psi}^{T}\mathbf{\Psi})^{-1}(\mathbf{\Psi}^{\text{T}}\mathbf{y}-\mathbf{q})$
(8)

Once $\boldsymbol{\alpha}$ is found, the denoised image can be estimated as

$\hat{\mathbf{x}}=\mathbf{\Psi}\boldsymbol{\alpha}$ (9)

A block diagram representation of the complete algorithm is shown in Fig. 2.
Note that in the figure we have assumed a two level UWT for simplicity. For
actual implementation we consider five levels of wavelet decomposition on
images of size $512\times 512$ which is row vectorized.

We have two fundamental operations in the GenRE algorithm: Wavelet analysis of
the image using undecimated Haar filter banks, which results in the matrix
$\mathbf{\Psi}$ and estimation of shrinkage parameter vector
$\boldsymbol{\alpha}$. Since this involves compute-intensive operations like
convolution, matrix inversion and multiplication, we should analyze them
separately and use every possible means to reduce computational complexity. A
complete analysis of this is done in the following sections.

##  IV Computational complexity and Memory usage of Haar Filter Bank

Two approaches to find the wavelet sub-bands are shown in Fig. 1 and Fig. 2.
They are in fact similar except that the filters in the later one are formed
by combining filters in a particular branch. For example, the first level Haar
filters are

$\displaystyle\mathbf{h}_{1}$
$\displaystyle=\frac{1}{2}\begin{bmatrix}1,1\end{bmatrix},\quad\mathbf{g}_{1}$
$\displaystyle=\frac{1}{2}\begin{bmatrix}-1,1\end{bmatrix},$
$\displaystyle\mathbf{\widetilde{h}}_{1}$
$\displaystyle=\frac{1}{2}\begin{bmatrix}1,1\end{bmatrix},\quad\mathbf{\widetilde{g}}_{1}$
$\displaystyle=\frac{1}{2}\begin{bmatrix}1,-1\end{bmatrix}.$ (10)

where $\mathbf{h}$ and $\mathbf{\widetilde{h}}$ are the refinement filters in
decomposition and recomposition phase, respectively. The corresponding wavelet
filters are $\mathbf{g}$ and $\mathbf{\widetilde{g}}$. Note that we have
normalized the filters by its length. Since we use UWT, filter for other
levels are found by up-sampling first level filters by $2^{j-1}$, where $j$
represents the level of filtering. Hence the filters for the second level are

$\displaystyle\mathbf{h}_{2}$
$\displaystyle=\frac{1}{2}\begin{bmatrix}1,0,1,0\end{bmatrix},\quad$
$\displaystyle\mathbf{g}_{2}=\frac{1}{2}\begin{bmatrix}-1,0,1,0\end{bmatrix},$
$\displaystyle\mathbf{\widetilde{h}}_{2}$
$\displaystyle=\frac{1}{2}\begin{bmatrix}1,0,1,0\end{bmatrix},\quad$
$\displaystyle\mathbf{\widetilde{g}}_{2}=\frac{1}{2}\begin{bmatrix}1,0,-1,0\end{bmatrix}.$
(11)

By using associativity of convolution and interchangeability of row and column
operations, we can combine the filters in each branch as

$\displaystyle\left.\begin{aligned}
\mathbf{h^{\prime}}_{1}&=\mathbf{h}_{1}=\frac{1}{2}\begin{bmatrix}1,1\end{bmatrix},\\\
\mathbf{h^{\prime}}_{2}&=\mathbf{h}_{1}*\mathbf{h}_{2}=\frac{1}{4}\begin{bmatrix}1,1,1,1\end{bmatrix},\\\
\mathbf{g^{\prime}}_{1}&=\mathbf{g}_{1}=\frac{1}{2}\begin{bmatrix}-1,1\end{bmatrix},\\\
\mathbf{g^{\prime}}_{2}&=\mathbf{h}_{1}*\mathbf{g}_{2}=\frac{1}{4}\begin{bmatrix}-1,-1,1,1\end{bmatrix},\\\
\mathbf{\widetilde{h^{\prime}}}_{1}&=\mathbf{\widetilde{h}}_{1}=\frac{1}{2}\begin{bmatrix}1,1\end{bmatrix},\\\
\mathbf{\widetilde{h^{\prime}}}_{2}&=\mathbf{\widetilde{h}}_{1}*\mathbf{\widetilde{h}}_{2}=\frac{1}{4}\begin{bmatrix}1,1,1,1\end{bmatrix},\\\
\mathbf{\widetilde{g^{\prime}}}_{1}&=\mathbf{\widetilde{g}}_{1}=\frac{1}{2}\begin{bmatrix}1,-1\end{bmatrix},\\\
\mathbf{\widetilde{g^{\prime}}}_{2}&=\mathbf{\widetilde{h}}_{1}*\mathbf{\widetilde{g}}_{2}=\frac{1}{4}\begin{bmatrix}1,1,-1,-1\end{bmatrix},\end{aligned}\right\\}$
(12)

where $*$ represents convolution. Filters for the other levels follow the same
structure but different lengths. The Haar filter bank in Fig. 2 uses these
combined filters.

So far, the filters we discussed are all in 1D domain. We can also do the
filtering in 2D domain, where the kernels are formed as outer products of 1D
filters. The corresponding individual kernels are

$\displaystyle\left.\begin{aligned}
\mathbf{HG}_{1}&=\frac{1}{4}\begin{bmatrix}-1&-1\\\
\phantom{-}1&\phantom{-}1\end{bmatrix},\\\
\mathbf{GH}_{1}&=\frac{1}{4}\begin{bmatrix}-1&\phantom{-}1\\\
-1&\phantom{-}1\end{bmatrix},\\\
\mathbf{GG}_{1}&=\frac{1}{4}\begin{bmatrix}\phantom{-}1&-1\\\
-1&\phantom{-}1\end{bmatrix},\\\
\mathbf{HH}_{1}&=\frac{1}{4}\begin{bmatrix}1&\phantom{-}1\\\
1&\phantom{-}1\end{bmatrix},\\\
\mathbf{HG}_{2}&=\frac{1}{4}\begin{bmatrix}-1&0&-1&0\\\
\phantom{-}0&0&\phantom{-}0&0\\\ \phantom{-}1&0&\phantom{-}1&0\\\
\phantom{-}0&0&\phantom{-}0&0\\\ \end{bmatrix},\\\
\mathbf{GH}_{2}&=\frac{1}{4}\begin{bmatrix}-1&0&1&0\\\ \phantom{-}0&0&0&0\\\
-1&0&1&0\\\ \phantom{-}0&0&0&0\\\ \end{bmatrix},\\\
\mathbf{GG}_{2}&=\frac{1}{4}\begin{bmatrix}\phantom{-}1&0&-1&0\\\
\phantom{-}0&0&\phantom{-}0&0\\\ -1&0&\phantom{-}1&0\\\
\phantom{-}0&0&\phantom{-}0&0\\\ \end{bmatrix},\\\
\mathbf{HH}_{2}&=\frac{1}{4}\begin{bmatrix}1&0&1&0\\\ 0&0&0&0\\\ 1&0&1&0\\\
0&0&0&0\\\ \end{bmatrix},\\\ \end{aligned}\right\\}$ (13)

and the combined kernels are

$\displaystyle\left.\begin{aligned}
\mathbf{HG}^{\prime}_{1}&=\mathbf{HG}_{1}=\frac{1}{4}\begin{bmatrix}-1&-1\\\
\phantom{-}1&\phantom{-}1\end{bmatrix},\\\
\mathbf{GH}^{\prime}_{1}&=\mathbf{GH}_{1}=\frac{1}{4}\begin{bmatrix}-1&\phantom{-}1\\\
-1&\phantom{-}1\end{bmatrix},\\\
\mathbf{GG}^{\prime}_{1}&=\mathbf{GG}_{1}=\frac{1}{4}\begin{bmatrix}\phantom{-}1&-1\\\
-1&\phantom{-}1\end{bmatrix},\\\
\mathbf{HH}^{\prime}_{1}&=\mathbf{HH}_{1}=\frac{1}{4}\begin{bmatrix}1&\phantom{-}1\\\
1&\phantom{-}1\end{bmatrix},\\\
\mathbf{HG}^{\prime}_{2}&=\mathbf{HH}_{1}*\mathbf{GG}_{2}=\frac{1}{16}\begin{bmatrix}-1&-1&-1&-1\\\
-1&-1&-1&-1\\\ \phantom{-}1&\phantom{-}1&\phantom{-}1&\phantom{-}1\\\
\phantom{-}1&\phantom{-}1&\phantom{-}1&\phantom{-}1\\\ \end{bmatrix},\\\
\mathbf{GH}^{\prime}_{2}&=\mathbf{HH}_{1}*\mathbf{GH}_{2}=\frac{1}{16}\begin{bmatrix}-1&-1&\phantom{-}1&\phantom{-}1\\\
-1&-1&\phantom{-}1&\phantom{-}1\\\ -1&-1&\phantom{-}1&\phantom{-}1\\\
-1&-1&\phantom{-}1&\phantom{-}1\\\ \end{bmatrix},\\\
\mathbf{GG}^{\prime}_{2}&=\mathbf{HH}_{1}*\mathbf{GG}_{2}=\frac{1}{16}\begin{bmatrix}\phantom{-}1&\phantom{-}1&-1&-1\\\
\phantom{-}1&\phantom{-}1&-1&-1\\\ -1&-1&\phantom{-}1&\phantom{-}1\\\
-1&-1&\phantom{-}1&\phantom{-}1\\\ \end{bmatrix},\\\
\mathbf{HH}^{\prime}_{2}&=\mathbf{HH}_{1}*\mathbf{HH}_{2}=\frac{1}{16}\begin{bmatrix}1&1&1&1\\\
1&1&1&1\\\ 1&1&1&1\\\ 1&1&1&1\\\ \end{bmatrix}.\end{aligned}\right\\}$ (14)

The corresponding recomposition kernels can be obtained by flipping the
kernels (LABEL:eq_Haar_coeff_2D_Decompose) and (14) both vertically and
horizontally. Kernels for the higher levels can be calculated in a similar way
and they hold the similar structure.

Now we have four sets of filters given by (10) – (14) and each of them can be
used to find UWT of the image independently. We need to select one set of
filters or a combination of them after analyzing computational complexity and
memory requirement. Across all the levels, individual filters have only 2 and
4 non-zero elements for 1D and 2D case, respectively. Here, we use a
convolution based filtering approach for hardware implementation. Furthermore,
the non-zero elements in both these filters are all unity. Hence we can
neglect the multiplication involved in the convolution. Combined filters, on
the other hand have all its elements as unity with a change of sign for half
of its elements. The change in the sign have some structure which is favorable
for recursive implementation and hence we can employ either convolution or
recursion based approaches. There is an additional multiplication with
normalization factor involved in both the cases which we neglect from
computation analysis, as it is a multiplication with power of 2 which can be
achieved by shift operation. A detailed analysis of the computations and
memory usage is done in the following sections.

###  IV-A Analysis of Individual Filters

The difference equations for 1D individual filters in decomposition can be
written as

$y(n)=x(n)+x(n-\frac{L}{2}),$ (15)

and

$y(n)=x(n-\frac{L}{2})-x(n),$ (16)

for refinement and wavelet filters, respectively,

Figure 3: Implementation of UWT using 1D individual filters. Here $x(n)$
represents row vectorized input image.

where L is the length of the filter in the particular filtering level. Here we
have one addition per output pixel. We need to apply the filter both on rows
and columns, as depicted in Fig. 3. Thus, there are a total of six additions
involved in the calculation of the four sub-bands in a level. The average
additions per output pixel is $6/4=1.5$. Further, the implementation of (15)
and (16) requires delays of length $\frac{L}{2}$. Since the image is row-
vectorized, we need $\frac{L}{2}\times 512$ sized buffers at the input of the
column filters to implement $\frac{L}{2}$ delay over columns. As a result, we
need one $\frac{L}{2}$-length buffer and two $\frac{L}{2}\times 512$ buffers
for the implementation of 1D filters in a particular level as shown in Fig. 3.

We can do a similar analysis in 2D. Using (LABEL:eq_Haar_coeff_2D_Decompose),
we write the equations for different sub-bands as

$\displaystyle\left.\begin{aligned}
\text{LH}&=&\text{D}-\text{A}+\text{C}-\text{B},\\\
\text{HL}&=&\text{D}-\text{A}-\text{C}+\text{B},\\\
\text{HH}&=&\text{D}+\text{A}-\text{C}-\text{B},\\\
\text{LL}&=&\text{D}+\text{A}+\text{C}+\text{B}.\end{aligned}\right\\}$ (17)

Here A, B, C and D represent
$\mathbf{I}(n_{1}-\frac{L}{2},n_{2}-\frac{L}{2}),\mathbf{I}(n_{1},n_{2}-\frac{L}{2}),\mathbf{I}(n_{1}-\frac{L}{2},n_{2})$
and $\mathbf{I}(n_{1},n_{2})$, respectively, where $\mathbf{I}$ is the input
image. Fig. 4 illustrates the implementation of these equations.

Figure 4: Implementation of UWT using 2D individual filters. Here A, B, C and
D represents
$\mathbf{I}(n_{1}-\frac{L}{2},n_{2}-\frac{L}{2}),\mathbf{I}(n_{1},n_{2}-\frac{L}{2}),\mathbf{I}(n_{1}-\frac{L}{2},n_{2})$
and $\mathbf{I}(n_{1},n_{2})$ respectively where $\mathbf{I}$ is the input
image and $x(n)$ is the row vectorized image.

Here we can share the intermediate results D$-$A, D$+$A, C$-$B and C$+$B.
Hence there are eight additions involved in the calculation of all the sub-
bands and the average additions per output pixel is $\frac{8}{4}=2$. Further,
we need only two buffers here, one $\frac{L}{2}$ and another
$\frac{L}{2}\times 512$ as compared to the 1D implementation. The bit-widths
of these buffers will also be less by one.

The filters in recomposition side are similar except for a flip. As a result,
the number of additions involved in filtering will also be the same. In
contrast to decomposition, inputs to recomposition filters are the different
sub-bands that we calculated in the decomposition phase. The results of these
filters are then combined to give a single low pass band from a filtering
level as shown in Fig. 1. Hence the average number of additions per output
pixel will increase to 6 in 1D and 12 in 2D implementations. The buffer usage
will also increase since we can not share the buffers as we did during
decomposition. Here, we may have to use four $\frac{L}{2}$ and
$\frac{L}{2}\times 512$ buffers for implementing the filters in a particular
level.

###  IV-B Analysis of Combined Filters

In contrast to the individual filters, combined filters have all elements with
unity magnitude, which results in $L-1$ and $L^{2}-1$ additions in 1D and 2D
domains, respectively. We can find the average additions per output pixel by
doing a similar analysis as in section IV-A with $L-1$ and $L^{2}-1$ replacing
their counter parts in the respective domains. This will result in
$\frac{6}{4}(L-1)$ and $L^{2}-1$ additions per output pixel respectively for
1D and 2D filters. As we have already mentioned, for these particular filters
we can use recursion in place of convolution to reduce calculations. For
instance, if we take the 1D wavelet filter in the fifth level, which is

$\mathbf{h^{\prime}}_{5}=\frac{1}{32}\begin{bmatrix}-1,-1,-1,\dots,1,1,1,\dots\end{bmatrix}_{1\times
32}$ (18)

we can implement it recursively as depicted in Fig. 5. Here, the initial value
is calculated using convolution and the remaining pixel values are found using
recursion. The recursive equations for this implementation are

$y(n)=y(n-1)-x(n-L)+x(n),$ (19)

and

$y(n)=y(n-1)-x(n-L)+2x(n-\frac{L}{2})-x(n),$ (20)

for refinement and wavelet filters, respectively, where L is the length of the
filter. There are three additions in the wavelet filter and two in the
corresponding refinement filter. Hence, calculation of LH and HL sub-bands
require five additions each. To find LL and HH sub-bands we need only five
more additions since we already have the partial results. The average
additions per output pixel are $(15)/4=3.75$. Implementation of recursive
filters for a single wavelet channel is shown in Fig. 7. Here, the requirement
of delays is similar to the individual filter case except that instead of
$\frac{L}{2}\times 512$ delays we need $L\times 512$ delays and we introduce a
new delay of size $1\times 512$ in the column filter side. Introduction of
this new delay is caused by the row-vectorization of the input image and
dependency of recursion on the previous result. Since the image is row-
vectorized, we apply the column filter row-wise, i.e., we calculate the first
result pixel in the first column and then the first pixel in the second column
and so on. Once the first row is finished we start with the second row. Now in
order to do a recursion, we need the result pixels in the previous row which
is stored in the $1\times 512$ buffer. Hence, altogether we need one
$L$-length, two $L\times 512$ and four $512$-length buffers for the
implementation of recursive filters in a particular level.

Now let’s examine the 2D filters. To illustrate recursion, we consider the
kernel for HH sub-band in level 3 which is given in Fig. 6(a).

Figure 5: Illustration of 1D recursive implementation of high pass filter in
fifth level of Haar filter bank. The initial value is found using convolution.

(a) 2D kernel for HH sub-band in level 3.

(b) Image pixels grouped corresponding to the filter structure.

Figure 6: Grouping of image pixels in level 3 to introduce 2D recursion.

Exploiting the structure of the filter, we can group the image pixels as shown
in Fig. 6(b), where A, B, C and D represent the sum of pixels in the group.
Now, HH sub-band can be calculated as

$\text{HH}=\text{D}+\text{A}-\text{C}-\text{B},$ (21)

which is same as (17). Here, the values A, B, C and D can be found by applying
a low-pass filter of size $\frac{L}{2}$x$\frac{L}{2}$ on the input. This low-
pass filter can be implemented by two 1D recursions given by (19) with length
of the filter as $\frac{L}{2}$ which results in four additions. Once the low
pass image is found, all the sub-bands can be calculated using (17). Hence,
the total number of additions required to find a sub-band is $4+3=7$. The
first four additions are common for all the sub-bands, hence the average
additions per output pixel is (4+3$\times$4)/4 = 4. Implementation of the
filter for a single wavelet channel is shown in Fig. 8. Here we need one
$\frac{L}{2}$ and two $\frac{L}{2}\times 512$ buffers which is less than the
1D recursive implementation. The bit-width requirements for these buffers are
also less.

Figure 7: Recursive implementation of filters in a single wavelet channel. The
recursion filter could be either high pass or low pass depending on the
channel.

Figure 8: Illustration of 2D recursive UWT for a single wavelet channel.

As already discussed in the previous section, the number of additions will
increase in the recombination side. For 1D implementation the average
additions will be 5 and for 2D, it will be 7.

We call the method of recursive implementation as recursive UWT (RUWT) for
further reference. Table I summarizes the average additions for filtering
level five. Since the delays are of large size, they are implemented using
Block RAMs (BRAM) available in FPGA. Table II gives the memory usage in terms
of number of 36K BRAMs in the fifth level of filtering. We have chosen level-
five for comparison since it has the maximum filter length across all the
levels, which is equal to 32.

TABLE I: Average additions per output pixel for filtering level five in different methods of implementing UWT Method | Decomposition | Recomposition
---|---|---
1D UWT | $1.50$ | $6.00$
2D UWT | $2.00$ | $12.00$
1D Convolution111Using combined filters | $46.50$ | $62.00$
2D Convolution11footnotemark: 1 | $1023.00$ | $1023.00$
1D RUWT | $3.75$ | $5.00$
2D RUWT | $4.00$ | $7.00$
TABLE II: Number of Block RAMs (36K) required for filtering level five in different methods of implementing UWT. Here we assume a bit-width of 16 for calculating the utilization. Method | Decomposition | Recomposition
---|---|---
1D UWT | 7.5 | 15
2D UWT | 4 | 15
1D RUWT | 17 | 32
2D RUWT | 15 | 30

We now select any one of the discussed methods or a combination of them for
the final implementation. UWT has the lowest number of additions among all.
However, we combine all the resulting sub-bands from a level to a single low
pass band during recombination. But to find the shrinkage factors, we need all
of them separately as shown in Fig. 2. Hence, UWT methods can not be used in
recombination side. This narrows down the list to UWT or RUWT for
decomposition and RUWT for recomposition. UWT has lesser number of
computations compared to RUWT. So we chose it for decomposition. We also have
to make a choice between 1D and 2D implementations. Since 2D methods have
lesser memory requirements we chose them. Finally, we select 2D UWT for
decomposition and 2D RUWT for recomposition.

###  IV-C Sub-band time alignment

We find the sixteen wavelet sub-bands which constitute the columns of the
matrix $\mathbf{\Psi}$ in (8) in parallel. The results will not be time
aligned since the kernel sizes are different.

Figure 9: When $\mathbf{GG}$ kernels in level 1 and 5 operate in parallel on
an image.

For example, if we consider the parallel application of filters in first and
fifth level on an image as shown in Fig. 9, the output of the decomposition
kernel in level one is available along with the second row, third element of
the input image, whereas the output from a filter in the fifth level will only
be available with the $16^{\text{th}}$ row, $17^{\text{th}}$ element. Thus
there is a delay of $16\times 512+16$ clock cycles between the corresponding
elements of sub-bands in the first and last level. Similar delays exist in
other levels too. However, to estimate the shrinkage factors in a pipelined
fashion, we need the sub-bands to be time aligned. This can be achieved by
delaying the inputs to the filters in levels one to four so that outputs from
all the levels are time aligned. This can also be done by delaying the outputs
which is not preferred since the output bit-width is more which leads to
increased memory utilization. However, there is a practical difficulty in
implementing this on the decomposition side as the input to a particular level
is the low pass output from the previous level.

Figure 10: Advance calculation of LL sub-bands in order to achieve time
alignment.

Now, if we delay the input to a particular level, inputs to the higher levels
will also get delayed by the same amount. Fig. 10 illustrates one way of
addressing this. Here we find the low pass sub-bands of levels one to four in
advance. These intermediate results are then properly delayed and applied to
the decomposition filters so that their outputs are time aligned.

###  IV-D Bit-width analysis and number format

We prefer fixed point implementation as it’s size, power consumption and
memory usage is less compared to floating point and can achieve more speed of
operation [14]. The input pixels of the image are eight bit unsigned numbers
and the corresponding Q format is Q8.0. Since the filtering involves equal
number of additions and subtractions, the range of the result from
decomposition filters after scaling will be $[-\frac{M}{2},\frac{M}{2}]$,
where $M=\max|I|$ and $I$ is the input pixel value. Hence the integer part of
the result can be represented in 7 bits and the extra bit becomes the sign
bit. The LL sub-band on the other hand have only additions, hence the output
after scaling will be unsigned with the same range as input and we need eight
bits to represent the integer part. The number of bits in the fractional part
is depended on the filter size and thus on the level of filtering. In the case
of recomposition filters, even though we have equal number of additions and
subtractions, the input to the filters are signed numbers. Here the integer
part of the input is represented by seven bits and hence same number of bits
are required to represent the integer part of the output as well. Table III
summarizes the change in Q format over different filtering levels except for
LL sub-band. For LL sub-bands the Q format is similar except that the integer
part is represented by 8 bits and also is unsigned.

TABLE III: Change of Q format for different sub-band outputs except for LL.
Input format is Q8.0(unsigned) and the outputs are signed.

Level | Scaling($\frac{1}{L})$ | Qm.n(Decomposition) | Qm.n(Recomposition)
---|---|---|---
1 | $2^{-2}$ | Q7.2 | Q7.4
2 | $2^{-4}$ | Q7.4 | Q7.8
3 | $2^{-6}$ | Q7.6 | Q7.12
4 | $2^{-8}$ | Q7.8 | Q7.16
5 | $~{}2^{-10}$ | Q7.10 | Q7.20

The bit-width at the output of the recombination filter varies from 12 bits at
the first level to 28 bits at the final level. We don’t need to keep all these
bits because during shrinkage we scale each of these sub-bands with
$\alpha_{i}$ which is always near unity or less and then combine and round to
the nearest eight bit integer number. Thus the fractional parts at farthest
positions will have less effect on the result. Hence we truncate the
fractional positions to eight bits as given in Table IV so that we have 16 bit
numbers at the maximum.

TABLE IV: Q format after truncation.

Level | Qm.n(Decomposition) | Qm.n(Recomposition)
---|---|---
1 | Q7.2 | Q7.4
2 | Q7.4 | Q7.6
3 | Q7.6 | Q7.6
4 | Q7.6 | Q7.6
5 | Q7.6 | Q7.6

##  V Estimation of Shrinkage Parameter

There are four matrix operations involved in the estimation of shrinkage
parameter vector $\alpha$ given in (8); a matrix multiplication, inversion of
resulting matrix, a matrix vector multiplication and a vector addition. Since
we use a five level wavelet decomposition, we will have sixteen wavelet sub-
bands which constitute the columns of the matrix $\mathbf{\Psi}$. We calculate
all these sixteen sub-bands in parallel and thus a row of this matrix becomes
available in each clock. We now write the matrix multiplication in (8) as an
outer product

$\mathbf{Q}=(\mathbf{\Psi}^{\text{T}}\mathbf{\Psi})=\mathbf{\Phi}\mathbf{\Phi}^{\text{T}}$
(22)

where $\mathbf{\Phi}=\mathbf{\Psi}^{\text{T}}$. Now the consecutive columns of
the matrix $\mathbf{\Phi}$ become available in successive clocks. The outer
product of a matrix can be found as summation of outer products of its columns
[15]. This outer product generation and accumulation can be realized using the
DSP cores available in FPGA. Since each of these columns have sixteen
elements, there are 256 parallel multiplications in the generation of outer
products. In this, 120 multiplications are redundant since $\mathbf{Q}$ is
symmetric. So we find the upper diagonal entries of $\mathbf{Q}$ first and
then reflect it along the diagonal to find the complete matrix. We need only
130 DSP cores to do this. We need another 16 DSP cores to do the matrix vector
multiplication in (8) and the vector addition can be implemented using LUTs.
We now rewrite (8) as

$\boldsymbol{\alpha}^{*}=\mathbf{Q}^{-1}\mathbf{c}$ (23)

where c = $(\mathbf{\Psi}^{\text{T}}\mathbf{y}-\mathbf{q})$. The vector
$\boldsymbol{\alpha}^{*}$ is the solution of the linear equation

$\mathbf{Q}\boldsymbol{\alpha}^{*}=\mathbf{c}$ (24)

Since $\mathbf{Q}$ is defined as $(\mathbf{\Psi}^{\text{T}}\mathbf{\Psi})$, it
is a positive semidefinite matrix. But in our denoising framework it can be
considered as a positive definite matrix as we assume the signal is corrupted
by noise which is uncorrelated. Hence we can use convex optimization
techniques to solve this equation instead of actually inverting the matrix to
find the solution, which require more resources. We use gradient descent
algorithm to solve (24) since it has lesser number of computations compared to
other optimization techniques [16]. Here the weight update equation can be
written as

$\boldsymbol{\alpha}_{i}=\boldsymbol{\alpha}_{i-1}+\mu(\mathbf{c}-\mathbf{Q}\boldsymbol{\alpha}_{i-1})$
(25)

The matrix vector multiplication in (25) can be implemented using 16 DSP
cores. Initial weight vector $\boldsymbol{\alpha}_{0}$ is taken as one, since
it is the ideal solution if there is no noise. The step size is chosen as
$\frac{1}{2^{13}}$ which is implemented using shift operation.

##  VI Hardware Implementation

The different blocks described in the previous section are now combined to
give the final architecture of the GenRE-Haar Image denoiser which is shown in
Fig. 11. Since we assume $512\times 512$ images, the matrix $\mathbf{\Psi}$
will have a size $262144\times 16$. This matrix is first used to find the
shrinkage $\boldsymbol{\alpha}$ and then to denoise the image. Hence we have
to store this until the shrinkage factor is calculated. Block RAMs can not be
used for this purpose since the matrix size is too large. To address this we
apply the input image twice. In the first iteration we calculate and store the
shrinkage factors. During the next iteration we denoise the image using these
stored shrinkage factors. The $\mathbf{Anlyse}/\overline{\mathbf{Denoise}}$
input serves this purpose. The performance of an FPGA design is determined by
its resource utilization and critical path delay which determines the
frequency of operation. Table V gives the resource utilization and maximum
operating frequency of the entire design. It has been observed that the
improvement in signal to noise ratio for full precision implementation is
negligible compared to 16 bit implementation.

Figure 11: Hardware architecture of GenRE-Haar Image denoiser. $\mathbf{x}$ and $\hat{\mathbf{x}}$ represents row vectorized input and output images respectively. TABLE V: Post implementation resource utilization of GenRE-Haar denoiser. The precision mentioned here refers to the precision of Haar filter bank. Resources | Full precision | 16 bit precision
---|---|---
#LUT | 30516 | 18397
#FF | 41349 | 23549
#BRAM36 | 150 | 132
#DSP | 168 | 168
FMax(MHz) | 169 | 183

##  VII Results

The proposed architecture is verified using standard test images corrupted
with noise and the performance of the design is evaluated using Peak Signal to
Noise Ratio (PSNR) which is a commonly used metric defined as

$\text{PSNR}=10\log_{10}(\frac{(\max(\mathbf{y}))^{2}}{\text{MSE}})\text{dB},$
(26)

$\text{where}\max(\mathbf{y})=\underset{i,j}{\max}~{}y(i,j)$. A comparison of
the PSNR obtained from MATLAB and FPGA implementations for images contaminated
with gaussian noise with $\sigma=25$ is shown in Fig. 12. Table VI – Table
VIII compares the performance for four different images and three different
noise distributions. We observe that the PSNR gain of FPGA implementation is
very close to that of MATLAB.

Another performance metric used is structural similarity (SSIM) index, which
measures the similarity between two images. It is defined as

$\text{SSIM}(I_{1},I_{2})=\frac{(2\mu_{I_{1}}\mu_{I_{2}}+c_{1})(2\sigma_{I_{1}I_{2}}+c_{2})}{(\mu_{I_{1}}^{2}+\mu_{I_{2}}^{2}+c_{1})(\sigma_{I_{1}}^{2}+\sigma_{I_{2}}^{2}+c_{2})},$
(27)

where
$\mu_{I_{1}}$ is the average of $I_{1}$,
$\mu_{I_{2}}$ is the average of $I_{2}$,
$\sigma_{I_{1}}$ is the variance of $I_{1}$,
$\sigma_{I_{2}}$ is the variance of $I_{2}$,
$\sigma_{I_{1}I_{2}}$ is the covariance of $I_{1}$ and $I_{2}$,
$k_{1}=0.01$ and $k_{2}=0.03$,
$L$ is the dynamic range of the pixel values, and
$c_{1}=(k_{1}L)^{2},\,\text{and}\,\,c_{2}=(k_{2}L)^{2}$ are two variables used
to stabilise the division with a weak denominator.

Table IX–XI compares the performance of the hardware in terms of SSIM index
and Table XII gives the execution time across different platforms. It has been
observed that the SSIM index we obtain from FPGA and MATLAB are very close and
the improvement in execution time is almost hundred times compared to the C
implementation.

(a) Clean Image.

(b) Noisy Image.
PSNR =
20.246 dB.

(c) MATLAB
reference
PSNR =
29.204 dB.
PSNR Gain =
8.958.

(d) FPGA result.
PSNR =
29.110 dB.
PSNR Gain =
8.864 dB.

(e) Clean Image.

(f) Noisy Image.
PSNR =
20.268 dB.

(g) MATLAB
reference.
PSNR =
27.621 dB.
PSNR Gain =
7.353 dB.

(h) FPGA result.
PSNR =
27.521 dB.
PSNR Gain =
7.253 dB

(i) Clean Image.

(j) Noisy Image.
PSNR =
20.236 dB.

(k) MATLAB
reference.
PSNR =
28.194 dB.
PSNR Gain =
7.958.

(l) FPGA result.
PSNR =
27.825 dB.
PSNR Gain =
7.589.

(m) Clean Image.

(n) Noisy Image.
PSNR =
20.171 dB.

(o) MATLAB
reference.
PSNR =
29.478 dB.
PSNR Gain =
9.307.

(p) FPGA result.
PSNR =
29.274 dB.
PSNR Gain =
9.104.

Figure 12: Comparison of PSNR improvement in MATLAB and FPGA implementations
for different test images contaminated with Gaussian noise with $\sigma=25$.
Here MATLAB implementation uses double precision and FPGA uses 16 bit
precision.

##  VIII Conclusion and future work

We presented an FPGA implementation of the GenRE based image denoising
algorithm in [13]. The proposed architecture gives a PSNR gain which is very
close to its MATLAB and C counter parts. The improvement in execution time is
around hundred times compared to the C implementation. We assumed images of
size $512\times 512$ for the current implementation and any higher dimensional
image can be processed as blocks of $512\times 512$. As a future work, we can
reduce the block size to $64\times 64$ which will reduce the resource
utilization and improve the performance. We can also introduce parallelism
between these $64\times 64$ blocks which will reduce the execution time.
Another improvement will be to extend the architecture for images contaminated
with multiplicative noise.

TABLE VI: Comparison of PSNR for Gaussian noise contamination with $\sigma=25$ Image | Input PSNR | PSNR (MATLAB) | PSNR (FPGA)
---|---|---|---
Lena | 20.246 | 29.204 | 29.110
Tank | 20.171 | 29.478 | 29.274
Boat | 20.268 | 27.621 | 27.512
House | 20.236 | 28.194 | 27.825
TABLE VII: Comparison of PSNR for Uniform noise contamination with $\sigma=25$ Image | Input PSNR | PSNR (MATLAB) | PSNR (FPGA)
---|---|---|---
Lena | 20.172 | 29.212 | 29.114
Tank | 20.168 | 29.453 | 29.264
Boat | 20.237 | 27.627 | 27.525
House | 20.244 | 28.212 | 27.848
TABLE VIII: Comparison of PSNR for Laplacian noise contamination with $\sigma=25$ Image | Input PSNR | PSNR (MATLAB) | PSNR (FPGA)
---|---|---|---
Lena | 20.352 | 29.186 | 29.074
Tank | 20.239 | 29.449 | 29.281
Boat | 20.377 | 27.631 | 27.532
House | 20.407 | 28.167 | 27.769
TABLE IX: Comparison of SSIM index for Gaussian noise contamination with $\sigma=25$ Image | Input SSIM | SSIM (MATLAB) | SSIM (FPGA)
---|---|---|---
Lena | 0.2740 | 0.7544 | 0.7474
Tank | 0.2618 | 0.6859 | 0.6746
Boat | 0.3481 | 0.7072 | 0.7034
House | 0.3247 | 0.7214 | 0.7107
TABLE X: Comparison of SSIM index for Uniform noise contamination with $\sigma=25$ Image | Input SSIM | SSIM (MATLAB) | SSIM (FPGA)
---|---|---|---
Lena | 0.2668 | 0.7530 | 0.7423
Tank | 0.2565 | 0.6866 | 0.6776
Boat | 0.3431 | 0.7060 | 0.7034
House | 0.3222 | 0.7213 | 0.7105
TABLE XI: Comparison of SSIM index for Laplacian noise contamination with $\sigma=25$ Image | Input SSIM | SSIM (MATLAB) | SSIM (FPGA)
---|---|---|---
Lena | 0.2842 | 0.7582 | 0.7495
Tank | 0.2730 | 0.6866 | 0.6791
Boat | 0.3614 | 0.7115 | 0.7075
House | 0.3403 | 0.7221 | 0.7063
TABLE XII: Comparison of execution times. System configuration and tools used for MATLAB and C implementation: Intel core i5-2300 @ 2.80GHz with 4 GB of RAM, Visual studio 2017 with Intel parallel studio for matrix operations Platform | Execution time (ms)
---|---
MATLAB | 1342
C | 380
FPGA | 3.6

## References

  * [1] Q. Corporation, “Understanding ccd read noise.” www.qsimaging.com/ccdnoise.html, 2008.  [Online; accessed 19-July-2017].
  * [2] M. H. Wilkinson and F. Schut, Digital image analysis of microbes: imaging, morphometry, fluorometry and motility techniques and applications.  John Wiley & Sons, 1998.
  * [3] G. E. Healey and R. Kondepudy, “Radiometric ccd camera calibration and noise estimation,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 16, no. 3, pp. 267–276, 1994.
  * [4] P. Gravel, G. Beaudoin, and J. A. De Guise, “A method for modeling noise in medical images,” IEEE Transactions on medical imaging, vol. 23, no. 10, pp. 1221–1232, 2004.
  * [5] V. S. Frost, J. A. Stiles, K. S. Shanmugan, and J. C. Holtzman, “A model for radar images and its application to adaptive digital filtering of multiplicative noise,” IEEE Transactions on pattern analysis and machine intelligence, no. 2, pp. 157–166, 1982.
  * [6] M. I. AlAli, K. M. Mhaidat, and I. A. Aljarrah, “Implementing image processing algorithms in fpga hardware,” in Applied Electrical Engineering and Computing Technologies (AEECT), 2013 IEEE Jordan Conference on, pp. 1–5, IEEE, 2013.
  * [7] A. Gabiger-Rose, M. Kube, R. Weigel, and R. Rose, “An fpga-based fully synchronized design of a bilateral filter for real-time image denoising,” IEEE Transactions on Industrial Electronics, vol. 61, no. 8, pp. 4093–4104, 2014.
  * [8] A. Gabiger, M. Kube, and R. Weigel, “A synchronous fpga design of a bilateral filter for image processing,” in Industrial Electronics, 2009. IECON’09. 35th Annual Conference of IEEE, pp. 1990–1995, IEEE, 2009.
  * [9] S. D. Dabhade, G. Rathna, and K. N. Chaudhury, “A reconfigurable and scalable fpga architecture for bilateral filtering,” IEEE Transactions on Industrial Electronics, vol. 65, no. 2, pp. 1459–1469, 2018.
  * [10] K. V. Hegde, V. Kulkarni, R. Harshavardhan, and S. S. David, “Adaptive reconfigurable architecture for image denoising,” in Parallel and Distributed Processing Symposium Workshop (IPDPSW), 2015 IEEE International, pp. 196–201, IEEE, 2015.
  * [11] S. Di Carlo, P. Prinetto, D. Rolfo, and P. Trotta, “Aidi: An adaptive image denoising fpga-based ip-core for real-time applications,” in Adaptive Hardware and Systems (AHS), 2013 NASA/ESA Conference on, pp. 99–106, IEEE, 2013\.
  * [12] M. Katona, A. Pižurica, N. Teslić, V. Kovačević, and W. Philips, “Fpga design and implementation of a wavelet-domain video denoising system,” in International Conference on Advanced Concepts for Intelligent Vision Systems, pp. 650–657, Springer, 2005.
  * [13] B. K. Subramanian, A. Gupta, and C. S. Seelamantula, “A distribution-independent risk estimator for image denoising,” in Proceedings of the Tenth Indian Conference on Computer Vision, Graphics and Image Processing, p. 52, ACM, 2016.
  * [14] A. Finnerty and H. Ratigner, “Reduce power and cost by converting from floating point to fixed point,” Tech. Rep. WP491 (v1.0), Xilinx, mar 2017.
  * [15] C. D. Meyer, Matrix analysis and applied linear algebra, vol. 71.  Siam, 2000.
  * [16] S. Bubeck et al., “Convex optimization: Algorithms and complexity,” Foundations and Trends® in Machine Learning, vol. 8, no. 3-4, pp. 231–357, 2015.