📄 255-261.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=255-261//-->
<!--UNASSIGNED1//-->
<!--UNASSIGNED2//-->
<CENTER>
<TABLE BORDER>
<TR>
<TD><A HREF="../ch05/253-254.html">Previous</A></TD>
<TD><A HREF="../ewtoc.html">Table of Contents</A></TD>
<TD><A HREF="261-266.html">Next</A></TD>
</TR>
</TABLE>
</CENTER>
<P><BR></P>
<H2><A NAME="Heading1"></A><FONT COLOR="#000077">Chapter 6<BR>Digital Audio Transform Recipes
</FONT></H2>
<P ALIGN="RIGHT">…then thus he cried,<BR>When I the process have in memory,<BR>How thou hast wearied me on every side…</P>
<P ALIGN="RIGHT"><I>—William Wordsworth, 1888</I></P>
<P>This chapter covers the finer points of one-dimensional digital signal processing. Chapter 5 showed some of the mathematical basis of the Fourier transform and pointed out some of its properties. This chapter explains how to implement the Fourier transform using uniformly sampled audio. We start with a description of the DFT (discrete Fourier transform) and IDFT (inverse discrete Fourier transform). We show how you can compute the DFT and IDFT in Java and demonstrate why such operations are usually slow. We discuss how to speed the DFT and IDFT using a class of algorithms known as the FFT (fast Fourier transform) and the IFFT (inverse fast Fourier transform).</P>
<P>We also demonstrate the computation of the PSD (power spectral density) mentioned in chapter 5. The remaining portion of the chapter is dedicated to applications of these transforms. The applications include filtering, windowing, pitch-shifting, and the spectral analysis of resampling.</P>
<H3><A NAME="Heading2"></A><FONT COLOR="#000077">The Discrete Fourier Transform</FONT></H3>
<P>Let
</P>
<P ALIGN="CENTER"><IMG SRC="images/06-01d.jpg"></P>
<P>(6.1)
</P>
<P>be the sampled version of the waveform <I>v(t)</I>, where <I>N</I> is the number of samples.</P>
<BLOCKQUOTE>
<P><FONT SIZE="-1"><HR><B>NOTE: </B>Equation 6.1 numbers from zero, rather than 1, to reflect the start point of arrays in Java.<HR></FONT>
</BLOCKQUOTE>
<P>Recall that the Fourier transform of <I>v(t)</I> was given by</P>
<P ALIGN="CENTER"><IMG SRC="images/06-02d.jpg"></P>
<P>(6.2)
</P>
<BLOCKQUOTE>
<P><FONT SIZE="-1"><HR><B>NOTE: </B>V(f) can exist only if <I>v</I>(<I>t</I>) is absolutely integrable:
<P ALIGN="CENTER"><B><I>j<SUB>n</SUB> j<SUB>n-1</SUB> ... j<SUB>1</SUB></I></B></P>
<P ALIGN="CENTER"><B><I>a<SUB>n</SUB> a<SUB>n-1</SUB> ... a<SUB>1</SUB></I></B></P>
<HR></FONT>
</BLOCKQUOTE>
<P ALIGN="CENTER"><IMG SRC="images/06-03d.jpg"></P>
<P>Recall also that the inverse Fourier transform of <I>V(f)</I> was given by</P>
<P ALIGN="CENTER"><IMG SRC="images/06-04d.jpg"></P>
<P>(6.3)
</P>
<P>To compute Equation 6.2 for the sampled waveform of Equation 6.1, we must perform DFT, which is given by</P>
<P ALIGN="CENTER"><IMG SRC="images/06-05d.jpg"></P>
<P>(6.4)
</P>
<P>Recall, from Chapter 5, Euler’s relation:</P>
<P ALIGN="CENTER"><IMG SRC="images/06-06d.jpg"></P>
<P>(6.4a)
</P>
<P>This permits the alternative formulation of the kernel of the transform in Equation 6.4b:</P>
<P ALIGN="CENTER"><IMG SRC="images/06-07d.jpg"></P>
<P>(6.4b)
</P>
<P>Recall that an even function is one that has the property <I>f</I> (-<I>x</I>) = <I>f</I> (<I>x</I>). Also, an odd function has the property <I>f</I> (-<I>x</I>) = -<I>f</I> (<I>x</I>). Because sin [theta] = -sin [theta] sine is an odd function. Also, because cos [theta] = -cos [theta] cosine is an even function. Using the odd-even function properties of sine and cosine, we rewrite Equation 6.4b:</P>
<P ALIGN="CENTER"><IMG SRC="images/06-08d.jpg"></P>
<P>(6.4c)
</P>
<P>Substituting Equation 6.4c into Equation 6.4, we obtain the form of the DFT used in the Java implementation:</P>
<P ALIGN="CENTER"><IMG SRC="images/06-09d.jpg"></P>
<P>(6.5)
</P>
<P>When implementing the Java program, we compute the real and imaginary parts of Equation 6.5 using the following:</P>
<!-- CODE SNIP //-->
<PRE>
for(int j = 0; j < N; j++) {
twoPijkOnN = twoPikOnN * j;
r_data[k] += v[j] * Math.cos( twoPijkOnN );
i_data[k] -= v[j] * Math.sin( twoPijkOnN );
}
r_data[k] /= N;
i_data[k] /= N;
</PRE>
<!-- END CODE SNIP //-->
<P>The complete routine for the DFT is given at the end of this section. Recall that the preceding loop is executed <I>N</I> times. The index, <I>k</I>, varies 0…<I>N</I>-1, because the DFT divides the spectrum into <I>N</I> parts. Each spectral division is accessed using the <I>frequency index</I>, <I>k</I>. The frequency, <I>f<SUB>k</SUB></I>, is computed from <I>k</I> using</P>
<P ALIGN="CENTER"><IMG SRC="images/06-10d.jpg"></P>
<P>(6.6)
</P>
<P>where the sampling period, [Delta], is given by</P>
<P ALIGN="CENTER"><IMG SRC="images/06-11d.jpg"></P>
<P>(6.7)
</P>
<P>An examination of Equations 6.5 and 6.6 shows that the spectrum is described by integral harmonics of <I>f<SUB>s</SUB> / N</I>. For example, suppose that the sampling rate, <I>f<SUB>s</SUB></I>, is 8,000 Hz and that the number of points is 2,048; then the smallest change in frequency that can be detected is given by</P>
<P ALIGN="CENTER"><IMG SRC="images/06-12d.jpg"></P>
<P>Therefore, integral harmonics of 4 Hz will be used to approximate <I>v(t)</I>.</P>
<P>The PSD (power spectral density) gives the power at a specific frequency index, <I>psd<SUB>k</SUB></I>. We compute the power at any <I>f<SUB>k</SUB></I> by summing the square of the real and imaginary components of the amplitude:</P>
<P ALIGN="CENTER"><IMG SRC="images/06-13d.jpg"></P>
<P>(6.8)
</P>
<P>The amplitude spectral density is given by the square root of the power spectral density.</P>
<BLOCKQUOTE>
<P><FONT SIZE="-1"><HR><B>CD-ROM: </B>An implementation of the DFT is shown next.<HR></FONT>
</BLOCKQUOTE>
<!-- CODE //-->
<PRE>
1. package lyon.audio;
2. import java.io.*;
3. import java.awt.*;
4. import grapher.Graph;
5. public class FFT extends Frame {
6. double r_data[] = null;
7. double i_data[] = null;
8. ...
public double[] dft(double v[]) {
int N = v.length;
double t_img, t_real;
double twoPikOnN;
double twoPijkOnN;
// how many bits do we need?
N=log2(N);
//Truncate input data to a power of two
// length = 2**(number of bits).
N = 1<<N;
double twoPiOnN = 2 * Math.PI / N;
// We truncate to a power of two so that
// we can compare execution times with the FFT.
// DFT generally does not need to truncate its input.
r_data = new double [N];
i_data = new double [N];
double psd[] = new double[N];
System.out.println(“Executing DFT on “+N+” points...”);
for(int k=0; k<N; k++) {
twoPikOnN = twoPiOnN *k;
for(int j = 0; j < N; j++) {
twoPijkOnN = twoPikOnN * j;
r_data[k] += v[j] * Math.cos( twoPijkOnN );
i_data[k] -= v[j] * Math.sin( twoPijkOnN );
}
r_data[k] /= N;
i_data[k] /= N;
psd[k] =
r_data[k] * r_data[k] +
i_data[k] * i_data[k];
}
return(psd);
}
</PRE>
<!-- END CODE //-->
<H4 ALIGN="LEFT"><A NAME="Heading3"></A><FONT COLOR="#000077">Bit Computations and a Log Review</FONT></H4>
<P>We have noticed that some students are confused by the use of logs in the following code:
</P>
<!-- CODE SNIP //-->
<PRE>
1. // how many bits do we need?
2. N=log2(N);
3. //Truncate input data to a power of two
4. // length = 2**(number of bits).
5. N = 1<<N;
</PRE>
<!-- END CODE SNIP //-->
<P>If lines 2 and 5 are obvious to you, you should probably skip to the next section. Otherwise, please read on! Logs are basic things, but people seem to forget them.
</P>
<P>The log function is defined by</P>
<P ALIGN="CENTER">if <I>x = a<SUP>Y</SUP></I> then <I>y</I> = log<SUB><I>a</SUB> x</I></P>
<P>For example, if <I>x</I> = 210, then 10 = log<SUB>2</SUB> <I>x</I>.</P>
<P>The following are known as the laws of logarithms:</P>
<P ALIGN="CENTER">log<SUB><I>a</I></SUB>(<I>xy</I>) = log<SUB><I>a</I></SUB><I>x</I> + log<SUB><I>a</I></SUB><I>y</I></P>
<P>For example, to find log<SUB>2</SUB> 4096 use the following:</P>
<P ALIGN="CENTER"><IMG SRC="images/06-14d.jpg"></P>
<P>Where ln is the natural log. Therefore,
</P>
<P ALIGN="CENTER"><IMG SRC="images/06-15d.jpg"></P>
<P>Also, to find log B<SUB><I>x</I></SUB> use this:</P>
<P ALIGN="CENTER"><IMG SRC="images/06-16d.jpg"></P>
<P>(6.9)
</P>
<P>Based on Equation 6.9, line 2 computes the number of bits needed, rounding down, via truncation, the integral number, <I>N</I>. To compute Equation 6.9 to the base 2, we use the following:</P>
<!-- CODE SNIP //-->
<PRE>
public static int log2(int n) {
return (int) (Math.log(n)/Math.log(2.0));
}
</PRE>
<!-- END CODE SNIP //-->
<P><BR></P>
<CENTER>
<TABLE BORDER>
<TR>
<TD><A HREF="../ch05/253-254.html">Previous</A></TD>
<TD><A HREF="../ewtoc.html">Table of Contents</A></TD>
<TD><A HREF="261-266.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 + -