%generate chirp signals
% wtitten for Jenélle Dowling by Bruce Land
% May 2010
%=== parameters to modify =====================================
% Sample rate
Fs = 44100 ; %standard audio sampling rate

%song
%total_song_duration = 1 ;% seconds
number_sweeps = 10 ;
sweep_time = 0.2 ;%seconds
sweep_gap_start = 0.05 ;%seconds
sweep_gap_end = .3 ;%seconds
sweep_method = 'quadratic' ; %Can be 'linear','quadratic', or 'logarithmic'
begin_freq = 1400 ;% Hz
end_freq = 4600 ;% Hz

% additive noise to simulate wind, etc
noise_gain = 0.1 ;% filtered white noise; set to zero for no noise
noise_filter = 500 ;% Hz -- noise lowpass frequency

% do you want a 'dominate freq'?
% 0 means no; 1 means yes
use_dom_freq = 1 ;
dom_freq = 3000 ;%Hz
dom_peak_size = 3 ;% height of spectral bump
dom_peak_width = .05 ;% width of spectral bump 0.01 to 0.25

% do you want a WAV file?
% 0 means no; 1 means yes
write_file = 1;

% do you want computation?
% 0 means no; 1 means yes
compute_psd = 1 ;
compute_spectrogram = 1;
max_analysis_freq = 6000 ;% HZ
spectrogram_window_size = 1024 ; % more means better f resolution
% less means more time resolution; 1024 is good balance
spectrogram_colors = 'gray'; %can be 'jet', 'hot' 'gray' 'bone' 
% 'autumn' 'spring' 'summer' 'copper' 'prism' 'cool'

% do you want to play song?
% 0 means no; 1 means yes
play_song = 1 ;

%=== end of parameters to modify ==============================

% basic timing
dt=1/Fs;
% compute the linear gap times
sweep_gaps = linspace(sweep_gap_start, sweep_gap_end, number_sweeps-1); 
% get the total length of the song
t=0:dt:(sweep_time*number_sweeps+sum(sweep_gaps)+.05) ;
song = zeros(size(t)) ;

% generate  chirps
chirp_time = 0:dt:sweep_time ;
song_segment = chirp(chirp_time, begin_freq, sweep_time, end_freq, sweep_method);

% surpress the beginning and ending click
click_window = ones(size(chirp_time));
click_window(1:400) = [1:400]/400 ;
click_window(end-399:end) = [400:-1:1]/400 ;
song_segment = song_segment .* click_window ;

% get the chirp start times/index
chirp_start_index(1) = 1;                
for i = 2:number_sweeps
    chirp_start_index(i) = fix((sweep_time+sweep_gaps(i-1))*Fs)...
        + chirp_start_index(i-1);
end

% now make the song
for i = 1:length(chirp_start_index)
    % song index list to be modified
    song_segment_index = ...
       [chirp_start_index(i):chirp_start_index(i)+length(chirp_time)-1];
    song(song_segment_index) = song_segment ;      
end

% add in noise (if any)
[b,a] = butter(2, noise_filter/(Fs/2));
song = song + noise_gain * filter(b,a,randn(size(t)));

%filter for dom freq
if use_dom_freq 
    [b,a] = butter(1,[dom_freq*(1-dom_peak_width), 
                dom_freq*(1+dom_peak_width)]/(Fs/2)) ;
    song = dom_peak_size*filter(b,a,song) + song  ;
end

%normalize song
song = song/max(song);

if play_song
    soundsc(song,Fs)
end

if compute_spectrogram
    figure(1);clf;
    spectrogram(song,[spectrogram_window_size],[],[],Fs, 'yaxis')
    set(gca, 'ylim', [0 max_analysis_freq])
    colormap(spectrogram_colors)
    title(['sweep song: time resolution=', ...
        num2str(spectrogram_window_size/Fs, '%1.3f'), 'sec'])
end

if compute_psd
    figure(2);clf;
    pwelch(song,[spectrogram_window_size],[],[],Fs)
    set(gca, 'xlim', [0 max_analysis_freq/1000])
    title('sweep song -- average power spectral density')
end

if write_file
    file_name = uiputfile('*.wav');
    wavwrite(song,Fs,16,file_name);
end