6 6.1 sound_wav_files flu00.wav.wav 44.1 khz 1/44100 spwave Text with Time spwave t T = t 44.1 khz t = 1 sec 44100 j t f j {f 0, f 1, f 2,, f 1 6.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) 5.1) ft) = 1 2 a 0 + a n cos 2nπt T + b n sin 2nπt ) T n=1 j t 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 + n=1 n=1 a n cos 2nπj 39 + b n sin 2nπj ) 6.1)
6.3 2n + p)πj p = 1, 2, cos 2n + p)πj = 2nπj + 2pπj = cos 2nπj 2n + p)πj, sin = sin 2nπj 6.1) n 6.1) f j = 1 2 a 0 + + + p=1 p=1 a p ) a p+1 cos 2πj + p=1 p=1 ) a p+2 cos 2 2πj + b p+1 ) p=1 b p+2 ) + sin 2πj sin 2 2πj + ) a p+ 1 cos p=1 2 1)πj ) + b p+ 1 sin p=1 2 1)πj 1 2 a 0 + a p a 0 a p+1 a p+2 p=1 a p+ 1 a 1 a 2 a 1 p=0 b p+1 b p+2 b p+ 1 b 1 b 2 b 1 p=0 p=0 f j = a 0 + 1 n=1 p=0 a n cos 2nπj + b n sin 2nπj ) p=0 j = 0, 1, 2,, 1) p=0 6.2) {f 0, f 1, f 2,, f 1 6.2) a 0 a 1 a 2 a 1 b 1 b 2 b 1 θ = 2π/ 1) i i = 1 p + 1 1 40
1 r = e ipθ S p = 0 S p 0 S 2) n j cos njθ cos n)jθ, sin njθ + sin n)jθ 3) 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 {f 0, f 1, f 2,, f 1 a n = 1 b n = 1 1 1 f k cos 2nπk f k sin 2nπk = 1 = 1 1 1 f k cos nkθ n = 0, 1, 2,, 1) f k sin nkθ n = 1, 2,, 1) θ = 2π 6.3) 6.2) f j 6.3) b 0 = 0 i i = 1 c n = a n ib n c n = 1 1 f k cos nkθ i sin nkθ) = 1 1 f k cos nkθ) + i sin nkθ)) = 1 6.2) a 0 + 1 n=1 a n cos njθ + b n sin njθ) = 1 c n e injθ = a n ib n )cos njθ + i sin njθ) 1 a n cos njθ + b n sin njθ) 6.3) f k e inkθ 6.4) = a n cos njθ + b n sin njθ) + ia n sin njθ b n cos njθ) 6.2) = 1 Re c n e injθ) = Re 1 c n e injθ ) 41
1 Re z) z c n e injθ 6.4) 1 c n e injθ = 1 1 1 f k e inkθ ) e injθ = 1 1 1 f k e ij k)nθ ) 1 e ij k)nθ 1 r = e ij k)θ + 1 j k 1 1) 1 c n e ijnθ = 1 1 1 f k e ij k)nθ ) = f j 6.2) f j 6.3) {f 0, f 1, f 2,, f 1 6.4 = 2M 6.2) 2) 3) n = 1, 2,, M 1 a n = a n, a n cos 2nπj = a n cos 2 n)πj b n = b n, sin 2Mπj f j = 2a 0) 2 + M 1 n=1 b n sin 2nπj = b n sin 2 n)πj = sin πj = 0 6.2) 2a n cos 2nπj + 2b n sin 2nπj ) + a M cos 2Mπj 6.5) j = 0, 1, 2,, 1) 2a n cos 2nπj + 2b n sin 2nπj ) 2nπj = 2s n cos ψ n s n = 6.5) f j = 2a 0) 2 a n2 + b n 2, + M 1 n=1 cos ψ n = a n s n, 2nπj 2s n cos ψ n sin ψ n = b n s n ) + a M cos 2Mπj 6.6) j = 0, 1, 2,, 1) 42
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 M T T a M { 2a 0, 2s 1, 2s 2,, 2s M 1, a M 6.7) 5000 4000 3000 2000 1000 0 0 1000 2000 3000 4000 5000 Figure 6.1: flute octave 4 Sol 5000 4000 3000 2000 1000 0 0 1000 2000 3000 4000 5000 Figure 6.2: harp octave 4 Sol Figure 6.1 flute octave 4 Sol 43
Figures 6.2 6.3 harp piano octave 4 Sol 5000 4000 3000 2000 1000 0 0 1000 2000 3000 4000 5000 Figure 6.3: piano octave 4 Sol 6.5 6.3) 6.3) 44.1 khz 0.5 = 22050 1 30 1965 2 = 2 k k Fast Fourier Transform FFT FFT FFT FFT ASA FFT 6.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 44
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] 6.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 6) es[0] es[1] es[2] es[m] 6.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) /* 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; 45
// 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 // 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; 46
// 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; 8) 1) flu04 flu04_fft.txt 1 T n = 0 0 0 T 2a 0 0 2a 0 2a 0 T 2s max s max s 1 s 2 s 1 2 n = M 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 fft_h_order.c fft_h_order./run.sh fft_h_order flu04.txt har04.txt pia04.txt tub04.txt vio04.txt ~/csc flute harp piano tuba violin octave 4 Do 47
Linux flute octave 4 Do fft_h_order flu04 flu04.txt flu04_fft.txt 2500 1800 1600 2000 1400 1200 1500 1000 1000 800 600 500 400 200 0 0 2000 4000 6000 8000 10000 12000 14000 0 0 2000 4000 6000 8000 10000 12000 14000 2000 1600 1800 1400 1600 1200 1400 1200 1000 1000 800 800 600 600 400 400 200 200 0 0 2000 4000 6000 8000 10000 12000 14000 0 0 2000 4000 6000 8000 10000 12000 14000 3000 2500 2000 1500 1000 500 0 0 2000 4000 6000 8000 10000 12000 14000 Figure 6.4: octave 4 Do flute harp piano tuba violin Figure 6.4 flute octave 4 Do gnuplot 300 15000 48
Linux harp piano tuba violin octave 4 Do gnuplot Figure 6.4 300 15000 6.6 {f 0, f 1, f 2,, f 1 6.2) a 0 a 1 a 2 a 1 b 1 b 2 b 1 6.3) a 0 a 1 a 2 a M b 1 b 2 b M 1 6.2) 6.5) 6.6) {f 0, f 1, f 2,, f 1 = 2M a n =a n b n = b n 6.6) s 1 s 2 s 1 s max 6.6) n = 1, 2,, M 1 f j = 2a 0) 2 + a M cos 2Mπj M 1 + n=1 2nπj 2s n cos ψ n ) j = 0, 1, 2,, 1) n = 1, 2,, M 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 m f [m] j f j idft.c { idft.c idft.c f [m] 0, f [m] 1,, f [m] 1 m 6.2) 6.5) { 6.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 49
2) a[0] a[1] a[2] a[m] b[1] b[2] b[m-1] 6.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); // 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++) 50
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 idft.c -lm -o idft = f [m] 0 m = 1, 2,, 9 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 51
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 6.7 Windows Q: csc.wav sound_wav_files sound_wav_files_cats_dogs sound_wav_files sound_wav_files_cats_dogs 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 4) Linux 5) fft_h_order FFT FFT 6) idft 7) 6) Windows 52
8) spwave Windows Linux & Windows 1) 8) 300 15000 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 set terminal postscript set output "filename.ps" replot set terminal postscript eps color set output "filename.eps" replot filename Windows Linux & Windows sound_wav_files_cats_dogs 1) 8) 53
200 10000 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 54