www.gusucode.com > matlab 案例源码 matlab代码程序 > matlab/SpectralAnalysisOfAWhaleCallExample.m

    %% Audio Signal
% This example analyzes the frequency spectrum of audio data.

%% 
% The file |bluewhale.au| contains audio data from a Pacific blue whale
% vocalization recorded by underwater microphones off the coast of
% California. The file is from the library of animal vocalizations
% maintained by the <http://www.birds.cornell.edu/brp/ Cornell University
% Bioacoustics Research Program>.

%%
% Because blue whale calls are so low, they are barely audible to humans.
% The time scale in the data is compressed by a factor of 10 to raise the
% pitch and make the call more clearly audible. Read and plot the audio
% data.  You can use the command |sound(x,fs)| to listen to the audio.
whaleFile = fullfile(matlabroot,'examples','matlab','bluewhale.au');
[x,fs] = audioread(whaleFile);

plot(x)
xlabel('Sample Number')
ylabel('Amplitude')

%%
% The first sound is a "trill" followed by three "moans". This example
% analyzes a single moan.  Specify new data that approximately consists of
% the first moan, and correct the time data to account for the factor-of-10
% speed-up.  Plot the truncated signal as a function of time.
moan = x(2.45e4:3.10e4);
t = 10*(0:1/fs:(length(moan)-1)/fs);

plot(t,moan)
xlabel('Time (seconds)')
ylabel('Amplitude')
xlim([0 t(end)])

%%
% The Fourier transform of the data identifies frequency components of the
% audio signal. In some applications that process large amounts of data
% with |fft|, it is common to resize the input so that the number of
% samples is a power of 2. This can make the transform computation
% significantly faster, particularly for sample sizes with large prime
% factors. Specify a new signal length |n| that is a power of 2, and use
% the |fft| function to compute the discrete Fourier transform of the
% signal. |fft| automatically pads the original data with zeros to increase
% the sample size.
m = length(moan);       % original sample length
n = pow2(nextpow2(m));  % transform length
y = fft(moan,n);        % DFT of signal

%%
% Adjust the frequency range due to the speed-up factor, and compute and
% plot the power spectrum of the signal.  The plot indicates that the moan
% consists of a fundamental frequency around 17 Hz and a sequence of
% harmonics, where the second harmonic is emphasized.
f = (0:n-1)*(fs/n)/10;
power = abs(y).^2/n;      

plot(f(1:floor(n/2)),power(1:floor(n/2)))
xlabel('Frequency')
ylabel('Power')