7 7.1 sound_wav_files flu00.wav.wav 44.1 khz 1/44100 spwave Text with Time spwave T > 0 t 44.1 khz t = 1 44100 j t f j {f 0, f 1, f 2,, f 1 = T t 7.2 T {f 0, f 1, f 2,, f 1 T ft) f j = fj t) j = 0, 1, 2,, 1; t = T ) ft) ft) 6.1) ft) = 1 2 a 0 + a n cos 2nπt T + b n sin 2nπt ) T j t j = 0, 1, 2,, 1 f j = fj t) = 1 2 a 0 + a n cos 2nπj t + b n sin 2nπj t ) T T = 1 2 a 0 + a n cos 2nπj 39 + b n sin 2nπj ) 7.1)
p = 1, 2, cos 2n + p)πj = cos 2nπj 2n + p)πj, sin = sin 2nπj 7.1) f j = 1 1 2 a ) 0 + a p + a n+p cos 2nπj p=1 p=0 1 + ) b n+p p=0 sin 2nπj 1 2 a 0 + a p a 0 a 1+p p=1 p=0 a 1+p a 1 a 1 b 1+p p=0 b 1+p b 1 b 1 p=0 p=0 f j = a 0 + 1 a n cos 2nπj + b n sin 2nπj ) j = 0, 1, 2,, 1) 7.2) 7.3 {f 0, f 1, f 2,, f 1 7.2) a 0 a 1 a 2 a 1 b 1 b 2 b 1 θ = 2π/ i i = 1 n 1 r = e inθ S 1) n = 0 S 2) n 0 S 3) j cos njθ cos n)jθ sin njθ + sin n)jθ 4) f 0 f 1 f 2 f 1 n a n = 1 1 f k cos nkθ, b n = 1 1 f k sin nkθ a n a n b n + b n 40
{f 0, f 1, f 2,, f 1 a n = 1 b n = 1 1 1 f k cos 2πnk f k sin 2πnk n = 0, 1, 2,, 1) n = 1, 2,, 1) 7.3) θ = 2π/ i i = 1 c n = a n ib n cos nkθ i sin nkθ = e inkθ 7.3) c n = 1 1 f k e inkθ n = 0, 1, 2,, 1) 7.4) b 0 = 0 7.2) a 0 + 1 = Re a n cos jnθ + b n sin jnθ) = 1 1 a n cos jnθ + b n sin jnθ) c n e ijnθ j = 0, 1, 2,, 1) 1 Re z z c n e ijnθ 7.4) 1 c n e ijnθ = 1 1 1 f k e inkθ ) e ijnθ = 1 1 1 f k e ij k)nθ ) 1 e ij k)nθ 1 r = e ij k)θ 1) 2) 1 c n e ijnθ = 1 1 1 f k e ij k)nθ ) = f j 7.3) {f 0, f 1, f 2,, f 1 = 2M 3) 4) n = 1, 2,, M 1 a n = a n, a n cos 2nπj = a n cos 41 2 n)πj
b n = b n, b n sin 2nπj = b n sin 2 n)πj sin 2Mπj f j = 2a 0) 2 + M 1 = sin πj = 0 7.2) 2a n cos 2nπj + 2b n sin 2nπj ) + a M cos 2Mπj 7.5) j = 0, 1, 2,, 1) 7.4 2a n cos 2nπj + 2b n sin 2nπj 2nπj = 2s n cos ψ n ) s n = a n2 + b n 2, cos ψ n = a n s n, sin ψ n = b n s n 7.5) f j = 2a 0) 2 + M 1 2nπj 2s n cos ψ n ) + a M cos 2Mπj 7.6) j = 0, 1, 2,, 1) 2nπj = 2π n T j t) 2s n {f 0, f 1, f 2,, f 1 n T 2a 0 1 T 2 T M 1 T M T a M { 2a 0, 2s 1, 2s 2,, 2s M 1, a M 7.7) Figure 2.1 flute octave 3 Do Figures 2.2 2.3 harp piano octave 3 Do 7.5 7.3) 7.3) 44.1 khz 0.5 42
= 22050 1 30 1965 2 = 2 k k Fast Fourier Transform FFT FFT FFT FFT ASA FFT 7.3) 1 30 FFT 2 3 FFT http://www.kurims.kyoto-u.ac.jp/~ooura/fft-j.html fft_h_order.c fft4g_h.c fft h order.c fft4g_h.c n = 1, 2,, M 1 1).txt flu04.txt flu04 2) T f 0 f f 0 f 3) FFT 2 < 2 k k 4) 2 k 4) 2 < 2 k f j = f 0 j = + 1, + 2,, 2 k 1) flu04 flu04_ext.txt 5) FFT rdft, 1, f) FFT rdft, 1, f) 1367 FFT rdft, 1, f) f[0] f[1] f[2] f[-1] f 0 f 1 f 2 f 1 rdft, 1, f) f[0] f[1] f[2] f[-1] 7.5) a 0 a 1 a 2 a M b 1 b 2 b 1 f[2n] = a n n = 0, 1, 2,, M 1 f[1] = a M f[2n+1] = b n n = 1, 2,, M 1 43
/* fft_h_order.c */ #include <math.h> #include <stdio.h> #include <string.h> #define MAX 32768) #define MAXH 16387) void rdftint, int, double *); //2^12 = 4096, 2^13 = 8192, 2^14 = 16384, 2^15 = 32768, 2^16 = 65536 int main) { int, tmp, M, i, j, k, n, flag, order[maxh+1]; double T, f[max+1], es[maxh+1]; FILE *fp; char file_name[100], file_in[100], file_out[100]; double zero = 0.0, tmp; // 1) printf"filename = "); scanf"%s", file_name); strcpyfile_in, file_name); strcatfile_in, ".txt"); fp = fopenfile_in, "r"); // 2) j = 0; flag = 1; while 1) { flag = fscanffp, "%lf %lf", &T, &f[j]); if flag == EOF) break; j = j+1; = j - 1; fclosefp); printf"origial: = %d T = %f\n",, T); f[0] = f[0] + f[])/2.0; f[] = f[0]; // 3) k = log)/log2.0); if pow2, k) < ) k = k + 1; tmp = pow2, k); // 4) 2 if tmp > ) { forj = + 1; j <= tmp; j++) f[j] = f[0]; T = T*tmp)/; = tmp; printf"added: = %d T = %f\n",, T); strcpyfile_out, file_name); strcatfile_out, "_ext.txt"); fp = fopenfile_out, "w"); forj = 0; j <= ; j++) fprintffp, "%f %f\n", j*t/, f[j]); fclosefp); // 5) FFT rdft, 1, f); // a[n] = f[2*n]*2.0/ n = 0, 1, 2,..., M-1) // a[m] = f[1]/ // b[n] = f[2*n+1]*2.0/ n = 1, 2, 3,..., M-1) // b[0] = b[m] = 0 44
// 6) M = /2; es[0] = fabsf[0])*2.0/; es[m] = fabsf[1])/; for n = 1; n < M; n++) { es[n] = sqrtf[2*n]*f[2*n]+f[2*n+1]*f[2*n+1])*2.0/; // 7) for n = 0; n <= M; n++) order[n] = n; for n = 2; n < M; n++) { for k = 1; k < n; k++) if es[k] < es[n]) break; tmp = es[n]; for i = n; i > k; i--) { es[i] = es[i-1]; order[i] = order[i-1]; es[k] = tmp; order[k] = n; // 8) strcpyfile_out, file_name); strcatfile_out, "_FFT.txt"); fp = fopenfile_out, "w"); n = 0; fprintffp, "%d %d %f %f %f %f %f %f\n",, n, n/t, f[n]*2.0/, zero, es[n], es[n]/es[1], T); fprintffp, "%d %d %f %f %f %f %f\n", M, M, M/T, f[1]/, zero, es[m], es[m]/es[1]); for k = 1; k < M; k++) { n = order[k]; fprintffp, "%d %d %f %f %f %f %f\n", k, n, n/t, f[2*n]*2.0/, f[2*n+1]*2.0/, es[k], es[k]/es[1]); return 0; 6) es[0] es[1] es[2] es[m] 7.7) a 0 s 1 s M 1 a M es[0] = 2a 0 es[n] = 2s n n = 1, 2,, M 1 es[m] = a M 7) n = 1, 2,, 1 8) 1) flu04 flu04_fft.txt 1 T n = 0 0 0 T 2a 0 0 2a 0 2a 0 2s max T s max s 1 s 2 s 1 2 n = M 45
M M M T a M 0 a M a M 2s max 3 n = 1, 2,, 1 n 2s 1 2s 2 2s M 1 k 2s n k+2 k n n T 2a n 2b n 2s n 2s n 3 n T 2s max 1 1 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 0 0 2000 4000 6000 8000 10000 12000 14000 0 2000 4000 6000 8000 10000 12000 14000 1 1 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 0 0 2000 4000 6000 8000 10000 12000 14000 0 2000 4000 6000 8000 10000 12000 14000 1 0.8 0.6 0.4 0.2 0 0 2000 4000 6000 8000 10000 12000 14000 Figure 7.1: octave 4 Do flute harp piano tuba violin csc10.tar Linux flu04.txt har04.txt pia04.txt 46
tub04.txt vio04.txt flute harp piano tuba violin octave 4 Do fft_h_order.c fft_h_order./run.sh Linux flute octave 4 Do fft_h_order flu04 flu04.txt flu04_fft.txt 2s max Figure 7.1 flute octave 4 Do gnuplot 300 15000 0.1 1.1 Linux harp piano tuba violin octave 4 Do gnuplot Figure 7.1 300 15000 0.1 1.1 {f 0, f 1, f 2,, f 1 7.2) a 0 a 1 a 2 a 1 b 1 b 2 b 1 7.3) a 0 a 1 a 2 a M b 1 b 2 b M 1 7.2) 7.5) 7.6) {f 0, f 1, f 2,, f 1 = 2M a n =a n b n = b n 7.6) n = 1, 2,, M 1 f j = 2a 0) 2 + a M cos 2Mπj M 1 + ) 2nπj 2s n cos ψ n j = 0, 1, 2,, 1) s max s 1 s 2 s 1 s n > s max /2 n f [1] j j = 0, 1, 2,, 1 f [1] j = 2a 0) + a M cos 2Mπj 2 + ) 2nπj 2s n cos ψ n. s n >s max /2 m = 2, 3, s n > s max /2 m n f [m] j j = 0, 1, 2,, 1 f [m] j = 2a 0) + a M cos 2Mπj 2 + ) 2nπj 2s n cos ψ n. s n >s max /2 m 47
{ idft.c idft.c f [m] 0, f [m] 1,, f [m] 1 m 7.2) { 7.5) 7.6) {f 0, f 1,, f 1 f [m] 0, f [m] 1,, f [m] 1 m = 1, 2,, 9 1) fft_h_order.c 8) flu04.txt flu04_fft 2) a[0] a[1] a[2] a[m] b[1] b[2] b[m-1] 7.5) a 0 a 1 a 2 a M b 1 b 2 b M 1 a[n] = 2a n n = 0, 1, 2,, M 1 a[m] = a M b[n] = 2b n n = 1, 2,, M 1 /* idfc.c */ #include <stdio.h> #include <math.h> #include <string.h> #define MAX 32768) #define MAXH 16387) int main ) { int, M, i, j, k, n, flag, order[maxh+1]; double T, f[max+1], a[maxh+1], b[maxh+1], es[maxh+1]; FILE *fp; char file_name[100], file_in[100], file_out[100], file_add[100]; int ktmp1, ktmp2; double dt, dtmp1, dtmp2, dtmp3, dtmp4, dtmp5; // 1) printf"filename = "); scanf"%s", file_name); strcpyfile_in, file_name); strcatfile_in, ".txt"); fp = fopenfile_in, "r"); // 2) fscanffp, "%d %d %lf %lf %lf %lf %lf %lf", &, &ktmp2, &dtmp1, &a[0], &dtmp2, &dtmp3, &dtmp4, &T); fscanffp, "%d %d %lf %lf %lf %lf %lf", &ktmp1, &ktmp2, &dtmp1, &dtmp2, &dtmp3, &dtmp4, &dtmp5); M = ktmp1; a[m] = dtmp2; k = 1; flag = 1; while 1) { flag = fscanffp, "%d %d %lf %lf %lf %lf %lf", &ktmp1, &n, &dtmp1, &dtmp2, &dtmp3, &dtmp4, &es[k]); if flag == EOF) break; order[k] = n; a[n] = dtmp2; b[n] = dtmp3; k = k+1; fclosefp); 48
// 3) dt = T/; for j = 0; j <= ; j++) f[j] = a[0]/2.0 + a[m]*cos2.0*m*m_pi*j/); k = 1; dtmp1 = 1.0; for i = 0; i <= 8; i++) { dtmp1 = dtmp1/2.0; printf"i = %d: ", i, k); while es[k]/es[1] >= dtmp1) { n = order[k]; for j = 0; j <= ; j++) f[j] = f[j] + a[n]*cos2.0*n*m_pi*j/) + b[n]*sin2.0*n*m_pi*j/); k = k + 1; if k > M) break; strcpyfile_out, file_name); sprintffile_add, "_IDFT%1d%05d.txt", i+1, k-1); strcatfile_out, file_add); fp = fopenfile_out, "w"); for j = 0; j <= ; j++) fprintffp, "%f %f\n", j*dt, f[j]); fclosefp); printf"%s includes until %d\n", file_out, k-1); return 0; 3) 1) flu04_fft s n > s max /2 n 9 s n > s max /2 2 n 21 s n > s max /2 3 n 40 s n > s max /2 9 n 827 { f [1] 0, f [1] 1,, f [1] flu04_fft_idft100009.txt { f [2] 0, f [2] 1,, f [2] flu04_fft_idft200021.txt { f [3] 0, f [3] 1,, f [3] flu04_fft_idft300040.txt { f [9] 0, f [9] 1,, f [9] flu04_fft_idft900827.txt f [m] idft.c idft cc -O2 idft.c -lm -o idft = f [m] 0 m = 1, 2,, 9 49
Linux & Windows flute octave 4 Do Linux idft flu04_fft flute octave 4 Do 9 flu04_fft_idft100009.txt flu04_fft_idft200021.txt flu04_fft_idft900827.txt 9 flu04.txt Windows spwave spwave flu04_fft_idft100009.txt Text with Time 44100 Linux & Windows Linux idft har04_fft harp octave 4 Do pia04_fft tub04_fft vio04_fft piano tuba violin octave 4 Do Windows spwave sound wav files cats dogs Windows sound_wav_files sound_wav_files_cats_dogs Windows.wav flu02.wav flute octave 3 Mi 9 1) 3) 8) Windows 5) 6) Linux 4) 7) Windows Linux 1) spwave flu02.wav 2) 0.743 3) Text with Time.txt 50
4) Linux 5) fft_h_order FFT FFT 6) idft 7) 6) Windows 8) spwave Windows Linux & Windows 1) 8) 300 15000 0.1 1.1 1 0 flu01.wav flu02.wav flu03.wav 3 1 1 flu05.wav flu06.wav flu07.wav 3 1 2 har01.wav har02.wav har03.wav 3 1 3 har05.wav har06.wav har07.wav 3 1 4 pia01.wav pia02.wav pia03.wav 3 1 5 pia05.wav pia06.wav pia07.wav 3 1 6 tub01.wav tub02.wav tub03.wav 3 1 7 tub05.wav tub06.wav tub07.wav 3 1 8 vio01.wav vio02.wav vio03.wav 3 1 9 vio05.wav vio06.wav vio07.wav 3 Windows Linux & Windows sound_wav_files_cats_dogs cat1.wav cat2.wav cat3.wav cat4.wav cat1.wav cat4.wav 3 1 1) 8) 200 10000 0.1 1.1 set terminal postscript set output "filename.ps" replot set terminal postscript eps color 51
set output "filename.eps" replot filename Windows Linux & Windows sound_wav_files_cats_dogs 1) 8) 200 10000 0.1 1.1 1 0 dog1.wav dog2.wav dog3.wav 3 1 1 dog4.wav dog b1.wav dog b2.wav 3 1 2 dog b3.wav dog b4.wav dog b5.wav 3 1 3 dog b6.wav dog b7.wav dog b8.wav 3 1 4 dog b9.wav dog bss1.wav dog bss2.wav 3 1 5 dog bss3.wav dog bss4.wav dog1.wav 3 1 6 dog2.wav dog3.wav dog4.wav 3 1 7 dog b1.wav dog b2.wav dog b3.wav 3 1 8 dog b4.wav dog b5.wav dog b6.wav 3 1 9 dog b7.wav dog b8.wav dog b9.wav 3 dog3.wav dog4.wav dog bss1.wav dog bss4.wav springer spaniel Windows Linux & Windows 1000 2000 52