📄 266-273.html
字号:
<HTML>
<HEAD>
<META name=vsisbn content="1558515682"><META name=vstitle content="Java Digital Signal Processing"><META name=vsauthor content="Douglas A. Lyon"><META name=vsimprint content="M&T Books"><META name=vspublisher content="IDG Books Worldwide, Inc."><META name=vspubdate content="11/01/97"><META name=vscategory content="Web and Software Development: Programming, Scripting, and Markup Languages: Java"><TITLE>Java Digital Signal Processing:Digital Audio Transform Recipes</TITLE>
<!-- HEADER --><STYLE type="text/css"> <!-- A:hover { color : Red; } --></STYLE><META NAME="ROBOTS" CONTENT="NOINDEX, NOFOLLOW">
<!--ISBN=1558515682//-->
<!--TITLE=Java Digital Signal Processing//-->
<!--AUTHOR=Douglas A. Lyon//-->
<!--PUBLISHER=IDG Books Worldwide, Inc.//-->
<!--IMPRINT=M & T Books//-->
<!--CHAPTER=6//-->
<!--PAGES=266-273//-->
<!--UNASSIGNED1//-->
<!--UNASSIGNED2//-->
<CENTER>
<TABLE BORDER>
<TR>
<TD><A HREF="261-266.html">Previous</A></TD>
<TD><A HREF="../ewtoc.html">Table of Contents</A></TD>
<TD><A HREF="273-278.html">Next</A></TD>
</TR>
</TABLE>
</CENTER>
<P><BR></P>
<H3><A NAME="Heading9"></A><FONT COLOR="#000077">Numeric Check of the DFT and IDFT</FONT></H3>
<P>A numeric check should be an integral part of every class. The <I>FFT</I> class contains a method called <I>testDFT</I> whose sole job is to verify the correctness of the DFT and IDFT implementation. With the number of samples set to 8, the testing method prints the input and output points for comparison.</P>
<!-- CODE //-->
<PRE>
public static void testDFT() {
int N = 8;
FFT f = new FFT(N);
double v[];
double x1[] = new double[N];
for (int j=0; j<N; j++)
x1[j] = j;
// take dft
f.dft(x1);
v = f.idft();
System.out.println(“j\tx1[j]\tre[j]\tim[j\t v[j]”);
for (int j=0; j < N; j++)
System.out.println(
j+”\t”+
x1[j]+”\t”+
f.r_data[j]+”\t”+
f.i_data[j]+”\t”+
v[j]);
}
</PRE>
<!-- END CODE //-->
<BLOCKQUOTE>
<P><FONT SIZE="-1"><HR><B>NOTE: </B>We print the intermediate DFT results to permit a detailed check against variations in implementation. Full data disclosure allows a baseline comparison among different implementations of the DFT. As we shall see in the following section, this data result is comforting, particularly when compared with the FFT data.<HR></FONT>
</BLOCKQUOTE>
<!-- CODE //-->
<PRE>
Executing IDFT on 8 points...
j x1[j] re[j] im[j] v[j]
0 0 3.5 0 -3.10862e-15
1 1 -0.5 1.20711 1
2 2 -0.5 0.5 2.00000
3 3 -0.5 0.207107 3
4 4 -0.5 0 4
5 5 -0.500000 -0.207107 5
6 6 -0.500000 -0.5 6
7 7 -0.5 -1.20711 7
</PRE>
<!-- END CODE //-->
<P>Although the input is not quite the same as the output, it is quite close.
</P>
<H3><A NAME="Heading10"></A><FONT COLOR="#000077">The FFT</FONT></H3>
<P>The FFT is a family of algorithms designed to speed the computation of the DFT:
</P>
<P ALIGN="CENTER"><IMG SRC="images/06-28d.jpg"></P>
<P>Direct computation of the DFT takes <I>O</I>(<I>N</I><SUP>2</SUP>) complex multiplications whereas the FFT takes <I>O</I>(<I>N</I> log <I>N</I>) complex multiplications. In this section we show an FFT algorithm known as the decimation in time, radix 2 FFT algorithm (also known as the Cooley-Tukey algorithm). The Cooley-Tukey algorithm is probably one of the most widely used of the FFT algorithms. The radix 2 property means that the number of samples to be transformed must be a power of two. The decimation in time means that the algorithm performs a recursive subdivision of the input sequence into its odd and even members. We are able to perform this subdivision as a result of the Danielson-Lanczos lemma:</P>
<P ALIGN="CENTER"><IMG SRC="images/06-29d.jpg"></P>
<P>Here is the proof of the Danielson-Lanczos lemma:
</P>
<P>Let</P>
<P ALIGN="CENTER"><IMG SRC="images/06-30d.jpg"></P>
<P>(6.13)
</P>
<P>so that</P>
<P ALIGN="CENTER"><IMG SRC="images/06-31d.jpg"></P>
<P>Substitute Equation 6.13 into Equation 6.4 to obtain
</P>
<P ALIGN="CENTER"><IMG SRC="images/06-32d.jpg"></P>
<P>(6.14)
</P>
<P>Separate Equation 6.14 into its odd and even components by altering the way the samples are indexed.</P>
<P ALIGN="CENTER"><IMG SRC="images/06-33d.jpg"></P>
<P>(6.15)
</P>
<P>Where Equation 6.15 shows summations operating over the odd and even indices. For example, if</P>
<P ALIGN="CENTER"><IMG SRC="images/06-34d.jpg"></P>
<P>then
</P>
<P ALIGN="CENTER"><IMG SRC="images/06-35d.jpg"></P>
<P>Factoring the exponents in Equation 6.15 yields
</P>
<P ALIGN="CENTER"><IMG SRC="images/06-36d.jpg"></P>
<P>The <I>W<SUP>k</SUP></I> term in the rightmost summation is not a function of the index, so</P>
<P ALIGN="CENTER"><IMG SRC="images/06-37d.jpg"></P>
<P>(6.16)
</P>
<P>To reflect the odd and even summations, Equation 6.16 is rewritten as</P>
<P ALIGN="CENTER"><IMG SRC="images/06-38d.jpg"></P>
<P>(6.17)
</P>
<P>Q.E.D.</P>
<P>The implications of Equation 6.17 are that we can divide the sequence into odd-and even-numbered samples. Thus, the Danielson-Lanczos lemma enables a divide and conquer algorithm to recursively split the sample sequence in half. The computational result of the Danielson-Lanczos lemma is that the <I>O</I>(<I>N</I><SUP>2</SUP>) DFT can be computed in <I>O</I>(<I>N</I> log <I>N</I>) time.</P>
<P>The Danielson-Lanczos lemma shows that a sequence must be divided into its odd and even subsets and these subsets must in turn be divided into their subsets. This process continues until we have only two members per subset. An illustration of this subdivision, for <I>N=8</I>, is shown in Figure 6.1.</P>
<P><A NAME="Fig1"></A><A HREF="javascript:displayWindow('images/06-01.jpg',264,162 )"><IMG SRC="images/06-01t.jpg"></A>
<BR><A HREF="javascript:displayWindow('images/06-01.jpg',264,162)"><FONT COLOR="#000077"><B>Figure 6.1</B></FONT></A> Decimation in time.</P>
<P>It is natural to implement the decimation in time using recursive calls with odd and even sets. It has been shown, however, that a recursive implementation is six times slower than a nonrecursive implementation [Gonzalez et al.]. To avoid a recursive approach to decimation, the bit-reversal algorithm (Cooley-Tukey) is used to perform the decimation in time (see Figure 6.2).
</P>
<P><A NAME="Fig2"></A><A HREF="javascript:displayWindow('images/06-02.jpg',357,197 )"><IMG SRC="images/06-02t.jpg"></A>
<BR><A HREF="javascript:displayWindow('images/06-02.jpg',357,197)"><FONT COLOR="#000077"><B>Figure 6.2</B></FONT></A> An example of how to decimate by bit reversal.</P>
<P>To arrive at the bit reversal, we implement a Java method in the <I>FFT</I> class:</P>
<!-- CODE SNIP //-->
<PRE>
int bitr(int j) {
int ans = 0;
for (int i = 0; i< nu; i++) {
ans = (ans <<1) + (j&1);
j = j>>1;
}
return ans;
}
</PRE>
<!-- END CODE SNIP //-->
<P>The <I>bitr</I> method works by linking two software shift registers, as shown in Figure 6.3.</P>
<P><A NAME="Fig3"></A><A HREF="javascript:displayWindow('images/06-03.jpg',191,61 )"><IMG SRC="images/06-03t.jpg"></A>
<BR><A HREF="javascript:displayWindow('images/06-03.jpg',191,61)"><FONT COLOR="#000077"><B>Figure 6.3</B></FONT></A> The <I>j</I> and <I>a</I> registers are linked with the + operator.
</P>
<P>After the decimation in time is performed, the balance of the computation is optimization hacks and housekeeping. For example, a simplification results from <I>Vk</I> being periodic in <I>N</I> so that <I>V<SUB>k+N</SUB></I> = <I>V<SUB>k</SUB></I>. Here is the proof:</P>
<P>Recall that the DFT is given by</P>
<P ALIGN="CENTER"><IMG SRC="images/06-39d.jpg"></P>
<P>so that
</P>
<P ALIGN="CENTER"><IMG SRC="images/06-40d.jpg"></P>
<P>Expanding the exponents and simplifying using
</P>
<P ALIGN="CENTER"><IMG SRC="images/06-41d.jpg"></P>
<P>with
</P>
<P ALIGN="CENTER"><IMG SRC="images/06-42d.jpg"></P>
<P>Thus,
</P>
<P ALIGN="CENTER"><IMG SRC="images/06-43d.jpg"></P>
<P>with
</P>
<P ALIGN="CENTER"><IMG SRC="images/06-44d.jpg"></P>
<P>(6.18)
</P>
<P>Q.E.D.</P>
<P>In addition, it can be shown that</P>
<P ALIGN="CENTER"><IMG SRC="images/06-45d.jpg"></P>
<P>(6.19)
</P>
<P>Here is the proof:</P>
<P>Recall Equation 6.13 is given by</P>
<P ALIGN="CENTER"><IMG SRC="images/06-46d.jpg"></P>
<P>So that
</P>
<P ALIGN="CENTER"><IMG SRC="images/06-47d.jpg"></P>
<P>with
</P>
<P ALIGN="CENTER"><IMG SRC="images/06-48d.jpg"></P>
<P>(6.20)
</P>
<P>This naturally leads to</P>
<P ALIGN="CENTER"><IMG SRC="images/06-49d.jpg"></P>
<P>Q.E.D.
</P>
<P>A further efficiency can be gained by the use of the recurrence relation:</P>
<P ALIGN="CENTER"><IMG SRC="images/06-50d.jpg"></P>
<P>(6.21)
</P>
<P>Here is the proof:</P>
<P ALIGN="CENTER"><IMG SRC="images/06-51d.jpg"></P>
<P>(6.22)
</P>
<P>Q.E.D.</P>
<P>The real part of Equation 6.22 is given by</P>
<P ALIGN="CENTER"><IMG SRC="images/06-52d.jpg"></P>
<P>so that
</P>
<P ALIGN="CENTER"><IMG SRC="images/06-53d.jpg"></P>
<P>(6.23)
</P>
<P>The imaginary part of Equation 6.22 is given by</P>
<P ALIGN="CENTER"><IMG SRC="images/06-54d.jpg"></P>
<P>so that
</P>
<P ALIGN="CENTER"><IMG SRC="images/06-55d.jpg"></P>
<P>(6.24)
</P>
<P>Equations 6.23 and 6.24 form the basis of the recurrence relationships, which enables the quick computation of the next <I>W<SUP>jk</SUP></I> based on the previous <I>Wjk</I>. An implementation of Equation 6.24 follows:</P>
<!-- CODE SNIP //-->
<PRE>
1. // (eq 6.23) and (eq 6.24)
2. wtemp = Wjk_r;
3. Wjk_r = Wj_r * Wjk_r - Wj_i * Wjk_i;
4. Wjk_i = Wj_r * Wjk_i + Wj_i * wtemp;
</PRE>
<!-- END CODE SNIP //-->
<P>Line 2 introduces <I>wtemp</I>, a temporary variable that enables the multiplication of the two complex numbers.</P><P><BR></P>
<CENTER>
<TABLE BORDER>
<TR>
<TD><A HREF="261-266.html">Previous</A></TD>
<TD><A HREF="../ewtoc.html">Table of Contents</A></TD>
<TD><A HREF="273-278.html">Next</A></TD>
</TR>
</TABLE>
</CENTER>
<hr width="90%" size="1" noshade><div align="center"><font face="Verdana,sans-serif" size="1">Copyright © <a href="/reference/idgbooks00001.html">IDG Books Worldwide, Inc.</a></font></div>
<!-- all of the reference materials (books) have the footer and subfoot reveresed --><!-- reference_subfoot = footer --><!-- reference_footer = subfoot --></BODY></HTML><!-- END FOOTER -->
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -