题面
众所周知的蜂鸣器
时间限制 : 10 Sec
空间限制 : 512 Mb
众所周知，Boshi 是个喜欢弹琴的妹子。
TA 经常在世界各地弹琴，比如 APIO/CTSC 赛场。
众所周知，[数据已删除] 是个喜欢唱歌的妹子。
TA 经常在世界各地唱歌，比如各种模拟考现场。
众所周知，Rayment 是个喜欢电音的妹子。
TA 经常在世界各地听电音，比如 Logic Gatekeeper。
众所周知，XZY 是个喜欢玩蜂鸣器的弱智。
TA 经常在世界各地玩蜂鸣器，然后被劝退。
最近 XZY 又在玩蜂鸣器了，TA 都玩了三年蜂鸣器了，从初中就开始玩了，玩到现在才玩出一个不着调的《极乐净土》，这让 TA 很难过。
TA 曾经想，如果能直接把下载下来的音乐直接转化成控制蜂鸣器的代码就好了。
控制蜂鸣器的代码长这样：
Beep(a, b);
表示蜂鸣器以 a 的频率振动（单位为赫兹），连续振动 b 毫秒。
TA 去世界各地查有没有这样的软件，但是众所周知，没有人会像 TA 这么无聊，所以没有。
TA 也很想自己写一个，但是因为自己太弱了所以写不出，TA 甚至都不知道音频文件的储存格式。
于是 TA 找啊找，找到了百度，找到了阿里巴巴，找到了四十大盗，找到了腾讯，找到了鲁迅，找到了枣树，找到了新中国。。。
最终 TA 在 Google 上找到了这样一段 c++代码可以读入.wav 格式的音乐（你可以不用读懂它在干什么的，等下 XZY 会向你说明的）：
#include <unistd.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <limits.h>
namespace WFR // Wave File Reader
{
double Length;
unsigned int num_samples;
int low_limit = 0l;
int high_limit = 0l;
int BUF[2][2000005];
int channels;
struct HEADER
{
unsigned char riff[4]; /* RIFF string */
unsigned int overall_size; /* overall size of file in bytes */
unsigned char wave[4]; /* WAVE string */
unsigned char fmt_chunk_marker[4]; /* fmt string with trailing null char */
unsigned int length_of_fmt; /* length of the format data */
unsigned int format_type; /* format type. 1-PCM, 3- IEEE float, 6 - 8bit A law, 7 - 8bit mu law */
unsigned int channels; /* no.of channels */
unsigned int sample_rate; /* sampling rate (blocks per second) */
unsigned int byterate; /* SampleRate * NumChannels * BitsPerSample/8 */
unsigned int block_align; /* NumChannels * BitsPerSample/8 */
unsigned int bits_per_sample; /* bits per sample, 8- 8bits, 16- 16 bits etc */
unsigned char data_chunk_header [4]; /* DATA string or FLLR string */
unsigned int data_size; /* NumSamples * NumChannels * BitsPerSample/8 - size of the next chunk that will be read */
};
unsigned char buffer4[4];
unsigned char buffer2[2];
struct HEADER header;
void Read()
{
int read = 0;
/* read header parts */
read = fread( header.riff, sizeof(header.riff), 1, stdin );
read = fread( buffer4, sizeof(buffer4), 1, stdin );
/* convert little endian to big endian 4 byte int */
header.overall_size = buffer4[0] | (buffer4[1] << 8) | (buffer4[2] << 16) | (buffer4[3] << 24);
read = fread( header.wave, sizeof(header.wave), 1, stdin );
read = fread( header.fmt_chunk_marker, sizeof(header.fmt_chunk_marker), 1, stdin );
read = fread( buffer4, sizeof(buffer4), 1, stdin );
/* convert little endian to big endian 4 byte integer */
header.length_of_fmt = buffer4[0] | (buffer4[1] << 8) | (buffer4[2] << 16) | (buffer4[3] << 24);
read = fread( buffer2, sizeof(buffer2), 1, stdin );
header.format_type = buffer2[0] | (buffer2[1] << 8);
read = fread( buffer2, sizeof(buffer2), 1, stdin );
channels = header.channels = buffer2[0] | (buffer2[1] << 8);
read = fread( buffer4, sizeof(buffer4), 1, stdin );
header.sample_rate = buffer4[0] | (buffer4[1] << 8) | (buffer4[2] << 16) | (buffer4[3] << 24);
read = fread( buffer4, sizeof(buffer4), 1, stdin );
header.byterate = buffer4[0] | (buffer4[1] << 8) | (buffer4[2] << 16) | (buffer4[3] << 24);
read = fread( buffer2, sizeof(buffer2), 1, stdin );
header.block_align = buffer2[0] | (buffer2[1] << 8);
read = fread( buffer2, sizeof(buffer2), 1, stdin );
header.bits_per_sample = buffer2[0] | (buffer2[1] << 8);
read = fread( header.data_chunk_header, sizeof(header.data_chunk_header), 1, stdin );
read = fread( buffer4, sizeof(buffer4), 1, stdin );
header.data_size = buffer4[0] | (buffer4[1] << 8) | (buffer4[2] << 16) | (buffer4[3] << 24);
/* calculate no.of samples */
num_samples = (8 * header.data_size) / (header.channels * header.bits_per_sample);
int size_of_each_sample = (header.channels * header.bits_per_sample) / 8;
/* calculate duration of file */
Length = (float) header.overall_size / header.byterate;
/* read each sample from data chunk if PCM */
if ( header.format_type == 1 ) /* PCM */
{
int i = 0, size_is_correct = 1;
char data_buffer[size_of_each_sample];
/* make sure that the bytes-per-sample is completely divisible by num.of channels */
int bytes_in_each_channel = (size_of_each_sample / header.channels);
if ( size_is_correct )
{
/* the valid amplitude range for values based on the bits per sample */
switch ( header.bits_per_sample )
{
case 8:
low_limit = CHAR_MIN;
high_limit = CHAR_MAX;
break;
case 16:
low_limit = SHRT_MIN;
high_limit = SHRT_MAX;
break;
case 32:
low_limit = INT_MIN;
high_limit = INT_MAX;
break;
}
for ( i = 1; i <= num_samples; i++ )
{
read = fread( data_buffer, sizeof(data_buffer), 1, stdin );
if ( read == 1 )
{
/* dump the data read */
unsigned int xchannels = 0;
int data_in_channel = 0;
for ( xchannels = 0; xchannels < header.channels; xchannels++ )
{
/* convert data from little endian to big endian based on bytes in each channel sample */
if ( bytes_in_each_channel == 4 )
{
data_in_channel = data_buffer[0] |
(data_buffer[1] << 8) |
(data_buffer[2] << 16) |
(data_buffer[3] << 24);
}
else if ( bytes_in_each_channel == 2 )
{
data_in_channel = data_buffer[0] |
(data_buffer[1] << 8);
}
else if ( bytes_in_each_channel == 1 )
{
data_in_channel = data_buffer[0];
}
BUF[xchannels][i] = data_in_channel;
/* check if value was in range */
if ( data_in_channel < low_limit || data_in_channel > high_limit )
printf( "**value out of range\n" );
}
}
else
{
printf( "Error reading file. %d bytes\n", read );
break;
}
} /* for (i =1; i <= num_samples; i++) { */
} /* if (size_is_correct) { */
} /* if (header.format_type == 1) { */
else puts( "Your Wave File isn't in PCM Format!" );
}
}
于是 TA 就能读入枣树，读入音乐了，TA 找到了人生的曙光，人民的希望，宇宙的奥秘，哲学的真理。
可是 TA 还是不会把音乐转化成 Beep，于是向你请教，并且告诉了你：
上面那段代码能这样调用：
WFR::Read()，无返回值，表示把.wav文件读入到WFR::BUF数组里WFR::Length，双精度小数型，表示读入的歌曲时长，单位为秒WFR::num_samples，无符号整型，表示整首歌的采样点个数WFR::low_limit，整形，表示振幅的下界（一般你是用不到的，这个跟采样位宽有关）WFR::high_limit，同上，表示振幅的上界（这里注意一下，振幅可以为负数，表示当前相位在初相位下方）WFR::BUF[i][j]，整形，表示声道 $i$的第 $j$个采样点的振幅，无单位，$i$从 $0$开始标号，$j$从 $1$开始标号WFR::channels，整形，表示声道个数上面的那段代码你可以自己随意修改。
如果你觉得上面那个
namespace太垃圾了自己写一个也是完全没有问题的。
当 XZY 讲完这些以后，爱弹琴的 Boshi 终于忍不住了，告诉了你：
实际上 XZY 的歌曲都很垃圾，最多只有两个声道，有的只有一个声道。
所以如果有两个声道的话你可以直接取
BUF[0][j]和BUF[1][j]的算术平均值作为第 $j$个采样点的振幅，这样就能把两个声道合成为一个声道了。而且如果你把一个声道的振幅数组画成柱状图，你会发现它其实就是声音信号的 “时域图”
对于单个音符，它的频率是单一的，固定的。
但是音乐里一般都有干扰，比如杂音。
XZY 应该只是想要你把杂音去掉，再把最明显的频率给 Beep 出来而已。
取出最明显的频率这种操作，估计要在 “频域图” 上动手脚了！
Boshi 讲完了这些，你恍然大悟，连连点头，欢天喜地，载兴载犇，僮仆欢迎，稚子候门，动手动脚，牛逼哄哄。
这时候，XZY 又说：
其实我的音乐都很简单的，杂音都很少，在底下的提示里我会写好每个歌曲的特点的！
于是 Boshi 把 XZY 的歌都听了一遍，发现不是很牛逼，于是现场演奏了一首 《两只老虎》，XZY 赶紧收藏了。
这时候一直躲在旁边的 [数据已删除] 坐不住了：[数据已删除] 天下第一！于是做现场的表演，演唱了一首 《？♂郝&@瀚 3#鸽》，在座的小朋友们都被 TA 精湛的演唱功底惊呆了，纷纷高兴地戳破了耳膜，XZY 也欢欣鼓舞一路火花加闪电地把美妙的歌声录了下来，赶紧收藏了。
屏幕前聪明的你，能否帮助聋了的 XZY 和大家完成这个任务呢？
输入格式
一个.wav 格式的文件，保证是 pcm 格式的。
另外：
此题需要开文件，从 wave.in 读入，输出到 wave.out。
请使用二进制方式打开文件，例如在 int main()的开始部分添加：
freopen("wave.in", "rb", stdin), freopen("wave.out", "w", stdout);
如果不用”rb” 方式打开文件可能会导致 Error reading file（本机也是！）！请注意！
输出格式
若干行，每行两个非负整数 $a$和 $b$，表示 Beep(a, b);
样例
请直接下载数据查看。
数据范围与提示
首先
输入文件是可以直接听的，下载下来以后把后缀名改为.wav，用音乐播放器打开就能听了。
你乐意用这个方法打表我很资磁
特殊问题
- 如果某一段时间原歌曲的声音很小，不足以判断出最明显的频率，请以前一段时间频率作为这一段时间的频率输出（若在歌曲开头处则频率直接设为 0）。
- 由于蜂鸣器都很菜，所以所有连续时长小于 $50$毫秒的 Beep 命令都无法发声（实际上如果运行的话会很奇怪，有的能发声有的不能），对于这种 Beep 命令 XZY 将视其 Beep 频率为 0（实际上就变成了 Sleep 了）
判断答案正误
XZY 写了个 SPJ 来判断你的输出是否正确
SPJ 会依次模拟执行你的 Beep 命令，如果出现以下情况则直接判 0 分：
- 命令个数超过 $10 ^ 5$
- 命令的总共执行时间（即播放时间）超过了 $6\times 10 ^ 5$毫秒
如果上面两条都没有发生，则 SPJ 会用该测试点的答案文件进行判分：
- 对于某一毫秒，如果两者的输出频率之差的绝对值小于 25（单位赫兹），或者标准答案这一毫秒是没有声音的（频率为 0），则判这一毫秒是正确的
- 对于整首歌曲，如果有 $40\%$及以上的时间是正确的（正确率高于 $40 \%$），则给满分
- 否则给出 “满分 $\times $正确率 $\div 40 \%$”，并且分数向下取整
数据范围
歌曲时长最大为 1 分 22 秒
采样点数不超过 $2 \times 10 ^ 6$个
其他基本不用操心
更详细的见下方 “歌曲特性”
歌曲特性
| 数据编号 | 歌曲名称 | 特性 |
|---|---|---|
| 1 | 两只老虎 | 由蜂鸣器演奏，XZY 用老人机录音 |
| 2 | 小酒窝 | 由蜂鸣器演奏，XZY 用老人机录音 |
| 3 | 极乐净土 | 由蜂鸣器演奏，XZY 用老人机录音 |
| 4 | 两只老虎 | Boshi 演奏，XZY 用老人机录音 |
| 5 | fc 坦克大战 BGM | 红白机风格，网易云下载，时长 4 秒 |
| 6 | 马里奥地上关 BGM | 红白机风格，XZY 用老人机录音 |
| 7 | 极乐净土 | 红白机风格，网易云下载 |
| 8 | 斗地主 | 红白机风格，网易云下载 |
| 9 | Only My Railgun | 红白机风格，网易云下载 |
| 10 | ？♂郝&@瀚 3#鸽 | Warning, nuclear missle launched! |
XZY 的老人机型号是 Nokia 216，单声道，采样率 8000（采样率就是每秒的采样点数目）
网易云音乐下载的都是双声道，采样率不会高于 11025
相关知识
第一组数据（两只老虎蜂鸣器版）的时域谱（横坐标对应时间，纵坐标对应振幅）：
第一组数据的频域谱（横坐标对应频率，纵坐标对应振幅，高峰表示突出的频率）：
题解
年度颓 (du) 废 (liu) 好题
感觉这个很好啊，应该可以算研究性学习了吧？
1. 采样与振幅
首先你要知道什么是 “采样”。
一首歌曲保存下来不可能是连续的，只可能是保存这个歌曲在某些时间点的信息。
一般情况下都是保存的 “振幅”
比如 “采样率” 是 8000，表示的是每秒钟的歌曲保存 8000 个采样点，每个采样点代表了 $\frac {1} {8000}$秒，每个采样点有一个有符号整数型值来储存它所代表的这一刻时间振膜的 “相位”
“相位” 说白了就是关于 “初始相位” 的位移，有正负号，“初始相位” 就是振膜不通电的时候的位置。
因为 $\frac {1} {8000}$秒十分短，所以这 $\frac {1} {8000}$秒内的相位可以视作都是相同的，因此可以用一个整数代替这 $\frac {1} {8000}$秒的相位。
显然，如果每秒采样次数越多，音质越高，也就是采样率越高音质越好，音乐文件也就越大。
（另外还有一个 “采样定理”，内容就是如果采样率高于原信号频率的两倍，那么这个采样就能完全代表原信号）
题目中读入的 WFR::BUF 就是每个采样点的整数信息，如果你用图形化界面输出这个数组的话就和题面给出的 “时谱图” 是一样的：
它的横坐标是时间（也就是采样点编号），纵坐标是振幅（也就是每个采样点的数值）。
可以看到中间有一条横线，那个是振幅为 0 的位置，也就是 “初始相位” 的位置。
在这条横线上面的表示振膜向上运动，下面的表示向下运动。
2. 时谱图——频谱图
现在我们有了时谱图的相关姿势，来讲下频谱相关姿势
什么是频谱？
比如我给你一个函数图像：
你可以把它看成一个时谱图对吧，横坐标是时间，纵坐标是振幅，振幅随着时间变化而周期性变化
我现在请你求出该函数的表达式，保证该函数由一堆正弦波叠加而成。
当然对正弦波熟悉的同学，肯定可以一下子看出这是 $f(x)= \sin (3 \pi x) + \sin (5 \pi x)$
但如果我给你的图像是：
这就比较无奈了对吧。。。（别试了这个有的频率都是小数了，而且我还在里面参了一个余弦波）
好了我们还是看回这个图像：
先说如果对这个图像进行频谱转换会得到什么：
你会得到这张图：
啥意思？
横坐标表示频率，纵坐标表示振幅。
它的意思是：你有一个振幅为 $1$的，频率为 $1.5$的周期函数（正弦波），还有一个振幅为 $1$，频率为 $2.5$的周期函数（正弦波）。
（其实正弦波和余弦波没啥区别，只不过是在时域上有一个相位差而已）
再看一个维基百科上的频谱分析的 3D 动画～
是不是很清晰明了了？
也就是说，你把你已知的时谱图转换成频谱图，就能分辨出原函数由哪些不同频率的正弦波/余弦波构成了！
当然还要说明一下，这个是对于周期函数说的，如果不是周期函数的话我们可以把它看成一个周期无限大的周期函数，然后再用正弦波/余弦波去拟合它即可。傅里叶大佬告诉我们，你用无限个正弦波余弦波总能拟合出来的（也比较显然对吧）
你看这个（正弦波叠加）：
现在的问题就是——怎么进行把时谱图转换成频谱图？
Warning! XZY 写了一整天也没想到怎么推出傅里叶级数公式，底下只是 XZY 写了一半，不忍心擦掉的部分。。。
好了，有请我们的装逼利器：欧拉公式！
$e ^ {ix} = \cos x + i \sin x$
$i$不必多讲吧，就是虚数单位，$i ^ 2 = -1$
所以。。。$e ^ {ix}$到底是什么呢？
我们把每一个三维坐标 $(\sin x, \cos x, x)$画到三维空间里～
下图中红轴表示 $\sin x$，绿轴表示 $\cos x$，蓝轴表示 $x$
（为了方便我把图形按红绿轴放大了 $5$倍，按蓝轴缩小了 $\pi$倍）
看！它升天了！它牛逼了！它载兴载犇，稚子候门，三径就荒，松菊犹存了！
这是什么意思呢？
红轴是 $\sin x$，也就是虚数轴
绿轴是 $\cos x$，也就是实数轴
蓝轴是 $x$，也就是时间轴
任意两个轴都构成了一个平面，我们先看看它在虚数轴和时间轴构成的平面上的投影（也就是在红轴和蓝轴构成的平面上的投影）：
哇，这不是一个有趣的、被旋转的正弦波吗？
再看看它在实数轴和时间轴构成的平面上的投影：
哇，这不是个有趣的、被旋转的余弦波吗？
（另外还有一个在虚实轴平面上的投影，显然就是个圆我都懒得放了）
嗯嗯，这可真是太棒了！
所以呢，你看啊：
$e ^ {ix} = \cos x + i \sin x$
$e ^ {-ix} = \cos x -i \sin x$
两式相加然后除个二得到：
$\cos x = \frac {e ^ {ix} + e ^ {-ix}} {2}$
欧耶，什么意思捏？
众所周知，$e ^ {ix}$表示一个随着时间变大而顺时针螺旋升天的奇妙卷卷。。。
那么 $e ^ {-ix}$表示的就是一个随着时间变大而逆时针螺旋升天的奇妙卷卷呀！
=￣ω￣=
$e ^ {-ix}$如下图所示（可以拿去和上面的 $e ^ {ix}$比♂较一下）：
那他们俩加起来是什么效果呢？
如下（红卷卷是 $e ^ {ix}$，蓝卷卷是 $e ^ {-ix}$，绿色的是 $2 \cos x$）
恩，就是这样的。
红色蓝色的虚部抵消了，就变成了绿色的波了。。。
（好吧其实上面这部分我只是想玩玩 Geogebra）
然后类似的可以证明：$\sin x = \frac {e ^ {ix} – e ^ {-ix}} {2i}$
好的，XZY 放弃了。。。XZY 太弱了
将上面两个关于 $\sin x$和 $\cos x$的式子带入到三角形式的傅里叶级数中就能得到复指数形式的傅里叶级数啦！！！！！！！！！！！！_(:з」∠)_
（又疯了一个）
。。。
XZY 真的不会了啊 (╬▔皿▔)
恩，反正式子长这样：
其中：
为什么？我怎么知道啊啊啊啊啊啊啊啊啊啊啊～！别问我啊！问维基百科啊！
反正你对时谱图数组进行一遍正向傅里叶变换就得到频谱图啦！！！
这时得到的频谱图是整首歌的频谱图，它的高峰代表的是整首歌中出现明显的频率，并不知道它们什么时候出现的，所以要把原本的一首歌分成很多很多小段（我是每 100 毫秒分一个小段），每个小段进行一遍傅里叶变换，得到这一段的最明显频率再输出。
其中有一些问题——傅里叶变换以后我们得到的到底是什么？
你对 $n$个采样点进行傅里叶变换，会得到 $n$个点的结果，有 $n$个复数。
其中第 $1$个复数（也就是数组的第 $0$位）代表的是 0Hz，第 $n$个复数（也就是数组的第 $n-1$位）代表的是 FsHZ，其中 Fs 表示你的采样率。
这很好理解，你傅里叶变换的数组的下标差是 1，而你采样率是 Fs，也就是数组里相邻的两项时间差是 $\frac {1} {Fs}$秒，因此实际上你是把时间轴扩大了 Fs 倍再进行时谱图转频谱图的，你只要最后把下标乘个 Fs 弄回来就行了。
然后关于某个频率的振幅，也就是频谱图的高度，应当虚数的什么属性代表，这个问题困扰了我很久。。。
实际上与虚数的模长有关。
如果原信号中的振幅（峰值）为 $P$，则 FFT 出来对应频率所对应的虚数的模长就是 $\frac {P\times N} {2}$，其中 $N$是 FFT 的长度。
大概是因为虚数模长代表了那个奇妙卷卷的半径。。。
（我记得网上还有个沙茶用实部作为振幅还发了一系列关于频谱分析的文章。。。）
最后还有一个问题：FFT 的时候一般都是取 $2 ^ k$格式的长度，那万一长度不整除怎么办呢？
实际上和多项式 FFT 没区别，在后面补零即可。经过试验，虽然补零会导致时域图变奇怪，但是补零不会太多，这样误差是完全可以忽略的。
这样就能愉快地打出代码啦！
分段进行 FFT，然后频谱分析一下，看哪个频率振幅最大输出。。。
复杂度基本上是 $O(n)$的
然后你会发现你的音乐有很多杂音（但其实因为 SPJ 为了能让大家过写的很水，所以应该可以直接 AC），这样你就需要优化了。。。
3. 优化
这里提供几个 XZY 的 XJB 优化方法：
- 时谱图中振幅低于最高振幅的 $0.1$倍的直接判为噪音，振幅设为 $0$，降低噪音
- 不考虑低于 $200$的频率，因为一般乐器不会低于 200Hz（反正数据中没有）
- 对于某一次 Beep，如果它的上一次 Beep 和下一次 Beep 频率相同且它的时长比较短，则判断当前这次 Beep 是杂音，把当前这次的频率设成上次的频率
- 最后合并一下相邻的且频率相似的 Beep
这样基本上就能比较美妙地唱出数据中的那 9 首歌和一首 [数据已删除] 演唱的锟斤拷了。
代码：
#include <unistd.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <limits.h>
namespace WFR // Wave File Reader
{
double Length;
unsigned int num_samples;
int low_limit = 0l;
int high_limit = 0l;
int BUF[2][2000005];
int channels;
struct HEADER
{
unsigned char riff[4]; /* RIFF string */
unsigned int overall_size; /* overall size of file in bytes */
unsigned char wave[4]; /* WAVE string */
unsigned char fmt_chunk_marker[4]; /* fmt string with trailing null char */
unsigned int length_of_fmt; /* length of the format data */
unsigned int format_type; /* format type. 1-PCM, 3- IEEE float, 6 - 8bit A law, 7 - 8bit mu law */
unsigned int channels; /* no.of channels */
unsigned int sample_rate; /* sampling rate (blocks per second) */
unsigned int byterate; /* SampleRate * NumChannels * BitsPerSample/8 */
unsigned int block_align; /* NumChannels * BitsPerSample/8 */
unsigned int bits_per_sample; /* bits per sample, 8- 8bits, 16- 16 bits etc */
unsigned char data_chunk_header [4]; /* DATA string or FLLR string */
unsigned int data_size; /* NumSamples * NumChannels * BitsPerSample/8 - size of the next chunk that will be read */
};
unsigned char buffer4[4];
unsigned char buffer2[2];
struct HEADER header;
void Read()
{
int read = 0;
/* read header parts */
read = fread( header.riff, sizeof(header.riff), 1, stdin );
read = fread( buffer4, sizeof(buffer4), 1, stdin );
/* convert little endian to big endian 4 byte int */
header.overall_size = buffer4[0] | (buffer4[1] << 8) | (buffer4[2] << 16) | (buffer4[3] << 24);
read = fread( header.wave, sizeof(header.wave), 1, stdin );
read = fread( header.fmt_chunk_marker, sizeof(header.fmt_chunk_marker), 1, stdin );
read = fread( buffer4, sizeof(buffer4), 1, stdin );
/* convert little endian to big endian 4 byte integer */
header.length_of_fmt = buffer4[0] | (buffer4[1] << 8) | (buffer4[2] << 16) | (buffer4[3] << 24);
read = fread( buffer2, sizeof(buffer2), 1, stdin );
header.format_type = buffer2[0] | (buffer2[1] << 8);
read = fread( buffer2, sizeof(buffer2), 1, stdin );
channels = header.channels = buffer2[0] | (buffer2[1] << 8);
read = fread( buffer4, sizeof(buffer4), 1, stdin );
header.sample_rate = buffer4[0] | (buffer4[1] << 8) | (buffer4[2] << 16) | (buffer4[3] << 24);
read = fread( buffer4, sizeof(buffer4), 1, stdin );
header.byterate = buffer4[0] | (buffer4[1] << 8) | (buffer4[2] << 16) | (buffer4[3] << 24);
read = fread( buffer2, sizeof(buffer2), 1, stdin );
header.block_align = buffer2[0] | (buffer2[1] << 8);
read = fread( buffer2, sizeof(buffer2), 1, stdin );
header.bits_per_sample = buffer2[0] | (buffer2[1] << 8);
read = fread( header.data_chunk_header, sizeof(header.data_chunk_header), 1, stdin );
read = fread( buffer4, sizeof(buffer4), 1, stdin );
header.data_size = buffer4[0] | (buffer4[1] << 8) | (buffer4[2] << 16) | (buffer4[3] << 24);
/* calculate no.of samples */
num_samples = (8 * header.data_size) / (header.channels * header.bits_per_sample);
int size_of_each_sample = (header.channels * header.bits_per_sample) / 8;
/* calculate duration of file */
Length = (float) header.overall_size / header.byterate;
/* read each sample from data chunk if PCM */
if ( header.format_type == 1 ) /* PCM */
{
int i = 0, size_is_correct = 1;
char data_buffer[size_of_each_sample];
/* make sure that the bytes-per-sample is completely divisible by num.of channels */
int bytes_in_each_channel = (size_of_each_sample / header.channels);
if ( size_is_correct )
{
/* the valid amplitude range for values based on the bits per sample */
switch ( header.bits_per_sample )
{
case 8:
low_limit = CHAR_MIN;
high_limit = CHAR_MAX;
break;
case 16:
low_limit = SHRT_MIN;
high_limit = SHRT_MAX;
break;
case 32:
low_limit = INT_MIN;
high_limit = INT_MAX;
break;
}
for ( i = 1; i <= num_samples; i++ )
{
read = fread( data_buffer, sizeof(data_buffer), 1, stdin );
if ( read == 1 )
{
/* dump the data read */
unsigned int xchannels = 0;
int data_in_channel = 0;
for ( xchannels = 0; xchannels < header.channels; xchannels++ )
{
/* convert data from little endian to big endian based on bytes in each channel sample */
if ( bytes_in_each_channel == 4 )
{
data_in_channel = data_buffer[0] |
(data_buffer[1] << 8) |
(data_buffer[2] << 16) |
(data_buffer[3] << 24);
}
else if ( bytes_in_each_channel == 2 )
{
data_in_channel = data_buffer[0] |
(data_buffer[1] << 8);
}
else if ( bytes_in_each_channel == 1 )
{
data_in_channel = data_buffer[0];
}
BUF[xchannels][i] = data_in_channel;
/* check if value was in range */
if ( data_in_channel < low_limit || data_in_channel > high_limit )
printf( "**value out of range\n" );
}
}
else
{
printf( "Error reading file. %d bytes\n", read );
break;
}
} /* for (i =1; i <= num_samples; i++) { */
} /* if (size_is_correct) { */
} /* if (header.format_type == 1) { */
else puts( "Your Wave File isn't in PCM Format!" );
}
}
#include <bits/stdc++.h>
#define NS (2000005)
#define PI (3.14159265358979323846)
using namespace std;
typedef complex<double> cpx;
int n, rev[NS * 3];
double data[NS], rate;
cpx d[NS * 3];
void FFT(cpx (&a)[NS * 3], int N, int bs, int f)
{
for (int i = 1; i < N; i += 1)
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (bs - 1));
for (int i = 0; i < N; i += 1)
if (i < rev[i]) swap(a[rev[i]], a[i]);
for (int l = 1; l < N; l <<= 1)
{
cpx w0(cos(PI / l), sin(PI / l) * f);
for (int i = 0; i < N; i += (l << 1))
{
cpx wn(1, 0), t1, t2;
for (int j = i; j < i + l; j += 1, wn *= w0)
t1 = a[j], t2 = wn * a[j + l], a[j] = t1 + t2, a[j + l] = t1 - t2;
}
}
if (f == -1) for (int i = 0; i < N; i += 1) a[i] /= N;
}
double mus[NS];
int cnt = 0;
int main(int argc, char const* argv[])
{
freopen("wave.in", "rb", stdin), freopen("wave.out", "w", stdout);
WFR::Read(), n = WFR::num_samples, rate = 1.0 * n / WFR::Length;
for (int i = 1; i <= n; i += 1)
{
for (int j = 0; j < WFR::channels; j += 1) data[i] += WFR::BUF[j][i];
data[i] /= (double)WFR::channels;
}
double mx_voc = *max_element(data + 1, data + 1 + n);
for (int i = 1; i <= n; i += 1)
if (data[i] < 0.1 * mx_voc)
data[i] = 0;
#define STEP 100
int sz = rate * STEP / 1000;
int N = 1, bs = 0; while (N <= sz) N <<= 1, bs++;
for (int i = 1; i + sz - 1 <= n; i += sz)
{
for (int j = 0; j < N; j += 1) d[j] = cpx(0, 0);
for (int j = i; j <= i + sz - 1; j += 1) d[j - i + 1].real(data[j]);
FFT(d, N, bs, 1);
double mx = 0, sum = 0, tot = 0; int mi = 0;
for (int j = 1; j <= (N >> 1); j += 1)
{
if (rate * j / N <= 200) continue;
sum += sqrt(norm(d[j])), tot++;
if (sqrt(norm(d[j])) > mx)
mx = sqrt(norm(d[j])), mi = j;
}
if (mx < sum / tot * 5) mus[++cnt] = 0;
else mus[++cnt] = rate * mi / N;
}
for (int i = 3; i <= cnt; i += 1)
if (abs(mus[i] - mus[i - 2]) <= 10 && mus[i] > 1e-4)
mus[i - 1] = mus[i];
for (int i = 1; i <= cnt;)
{
double F = mus[i]; int sum = STEP; i++;
while (i <= cnt && (abs(mus[i] - F) <= 10 || mus[i] < 1e-4))
{
if (mus[i] < 1e-4) mus[i] = F;
sum += STEP, F = (mus[i] + F) * 0.5, i++;
}
printf("%d %d\n", (int)F, sum);
}
return 0;
}
