p = 1, 2, cos 2n + p)πj = cos 2nπj 2n + p)πj, sin = sin 2nπj 7.1) f j = 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 +

Size: px

Start display at page:

Download "p = 1, 2, cos 2n + p)πj = cos 2nπj 2n + p)πj, sin = sin 2nπj 7.1) f j = 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 +"

しほこかりこめ
7 years ago
Views:

1 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 = 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)

2 p = 1, 2, cos 2n + p)πj = cos 2nπj 2n + p)πj, sin = sin 2nπj 7.1) f j = 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 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

3 {f 0, f 1, f 2,, f 1 a n = 1 b n = 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 = 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θ = f k e inkθ ) e ijnθ = f k e ij k)nθ ) 1 e ij k)nθ 1 r = e ij k)θ 1) 2) 1 c n e ijnθ = 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

4 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 harp piano octave 3 Do ) 7.3) 44.1 khz

5 = = 2 k k Fast Fourier Transform FFT FFT FFT FFT ASA FFT 7.3) 1 30 FFT 2 3 FFT 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

6 /* 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 = 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

7 // 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 = T 2a 0 0 2a 0 2a 0 2s max T s max s 1 s 2 s 1 2 n = M 45

8 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 Figure 7.1: octave 4 Do flute harp piano tuba violin csc10.tar Linux flu04.txt har04.txt pia04.txt 46

9 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 Linux harp piano tuba violin octave 4 Do gnuplot Figure {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

10 { 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

11 // 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_idft txt { f [2] 0, f [2] 1,, f [2] flu04_fft_idft txt { f [3] 0, f [3] 1,, f [3] flu04_fft_idft txt { f [9] 0, f [9] 1,, f [9] flu04_fft_idft txt f [m] idft.c idft cc -O2 idft.c -lm -o idft = f [m] 0 m = 1, 2,, 9 49

12 Linux & Windows flute octave 4 Do Linux idft flu04_fft flute octave 4 Do 9 flu04_fft_idft txt flu04_fft_idft txt flu04_fft_idft txt 9 flu04.txt Windows spwave spwave flu04_fft_idft txt Text with Time 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) ) Text with Time.txt 50

13 4) Linux 5) fft_h_order FFT FFT 6) idft 7) 6) Windows 8) spwave Windows Linux & Windows 1) 8) flu01.wav flu02.wav flu03.wav flu05.wav flu06.wav flu07.wav har01.wav har02.wav har03.wav har05.wav har06.wav har07.wav pia01.wav pia02.wav pia03.wav pia05.wav pia06.wav pia07.wav tub01.wav tub02.wav tub03.wav tub05.wav tub06.wav tub07.wav vio01.wav vio02.wav vio03.wav 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) set terminal postscript set output "filename.ps" replot set terminal postscript eps color 51

14 set output "filename.eps" replot filename Windows Linux & Windows sound_wav_files_cats_dogs 1) 8) dog1.wav dog2.wav dog3.wav dog4.wav dog b1.wav dog b2.wav dog b3.wav dog b4.wav dog b5.wav dog b6.wav dog b7.wav dog b8.wav dog b9.wav dog bss1.wav dog bss2.wav dog bss3.wav dog bss4.wav dog1.wav dog2.wav dog3.wav dog4.wav dog b1.wav dog b2.wav dog b3.wav dog b4.wav dog b5.wav dog b6.wav 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

6 6.1 sound_wav_files flu00.wav.wav 44.1 khz 1/44100 spwave Text with Time spwave t T = N t N 44.1 khz t = 1 sec j t f j {f 0, f 1, f 2,, f N 1

$6 6.1 sound_wav_files flu00.wav.wav 44.1 khz 1/44100 spwave Text with Time spwave t T = N t N 44.1 khz t = 1 sec j t f j {f 0, f 1, f 2,, f N 1$ 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,,