Преброяване на цвърченията на Neoconocephalus ensiger

Зареждане на аудио информацията

$HistoryLength=0;

sample=;

snd=sample[[1]][[1]][[1]];

srate=44100;

duration=N@((Length@snd)/srate)

Честотен анализ

Да видим звука. В следната графика по y е силата, по x – времето (в семпли).

ListPlot@snd

Използвайки бързата трансформация на Фурие (много) бързо ще получим честотната картина за звука.

sndFR=FourierDCT@snd;

ListPlot@sndFR

По y отново е силата, но по x вече е честотата (в Hz). Въпреки това, числата не изглеждат точни. Трябва да рескалираме абсцисата.

По теория, долната граница се определя от:

lowfreq=(1/duration)

А горната:

highfreq=(srate/2)

Тогава стъпката е:

h=((highfreq-lowfreq)/(Length@sndFR))

Две функции за преминаване от брой семпли в честоти и обратно:

sample2freq[i_]:=Round[i h]

freq2sample[f_]:=Round@(f/h)

Сега вече можем да получим коректната картина:

ListPlot@
Table[
{lowfreq+sample2freq@i,sndFR[[i+1]]}
,{i,0,Length@sndFR-1}
]

Четейки изследванията на Pierce, 1948, както и други такива, очакванията са честотите да са главно в региона 12–16 KHz. При нас това не е така понеже записът включва и други звуци (жаби и др.). Нека филтрираме излишните звуци

Честотен филтър

fltfreq[data_,p_,q_]:=
Module[
{
res=data,
myp=p,
myq=q
},
myp=freq2sample@p;
myq=freq2sample@q;
If[myp>1,
res[[1;;myp]]=0
];
If[myq<Length@data,
res[[myq;;]]=0
];
res
]

Да чуем изчиствайки звуците извън региона 12–16 KHz.

Sound@SampledSoundList[FourierDCT@fltfreq[sndFR,12000,16000],srate]

Нека сме по-точни. По Pierce знаем, че очакваните честоти са 13 700 Hz.

Sound@SampledSoundList[FourierDCT@fltfreq[sndFR,13690,13710],srate]

Това е, което търсим.

snd2=FourierDCT@fltfreq[sndFR,13690.0,13710];

Нека погледнем отблизо началото. Виждаме насложените вълни.

ListPlot[snd2[[;;4000]]]

Нека получим само една вълна.

Слепващ филтър

fltjoin[data_,step_]:=
Module[
{
res,
n=Floor@((Length@data)/step),
i
},
res=Array[1&,n+1];
For[i=0,i<n,i++,
res[[i+1]]=Max@data[[i step+1;;(i+1)step ]];
];
res[[n+1;;]]=Max@data[[n step+1;;]];
res
]

snd3=fltjoin[snd2,100];

Да разгледаме слепената вълна отблизо:

Manipulate[
ListPlot[snd3[[l;;l+40]]
,PlotRange->{{0,40},{0,0.0065}}]
,{l,1,Length@snd3-40,1}
]

Отчитане и броене на пиковете

Сега остава да отчетем пиковете.

Чрез шаблон (бавно, но директно)

Чрез шаблон това можем да направим най-директно, но пък е изключително не ефективно.

Timing[snd4=snd3[[;;140]]//.{m___,n__,i_,p__,q___}/;n<i&&p<i&&i≠2->{m,n,2,p,q};
]

Забелязва се как върховете са били "изкоренени".

ListPlot@snd4

Колко пика има в този интервал?

Cases[snd4,2]

Чрез цикъл

Чрез цикъл нещата вече са мигновени.

fltpeaks[data_]:=
Module[
{
peak,
count=0,
state=0,
i
},
For[i=1,i<Length@data,i++,
If[data[[i]]<data[[i+1]],state=1];
If[data[[i]]>data[[i+1]],
If[state==1,
(* т.е. имаме връх *)
peak[++count]={i,data[[i]]};
];
state=-1;
];
];
Array[peak,count]
]

Да видим колко, къде и какви върхове има.

peaks=fltpeaks@snd3

Да изобразим получения резултат.

sndpeaks=Array[0&,Length@snd3];

For[i=1,i≤Length@peaks,i++,
sndpeaks[[peaks[[i,1]]]]=peaks[[i,2]];
]

ListPlot[{snd3,sndpeaks},Joined->True]

Броят цвърчения е:

((Length@peaks)/duration)