audio_spectrum.m

filenames = {'samples/teacher/super.wav','samples/teacher/yume_deem.wav','samples/teacher/yume2005_deem.wav'};
%'samples/panic.wav','samples/super.wav', 'samples/yume.wav','samples/yume2005.wav','samples/single.mp3'
%filenames = {'old/good_rg.wav','old/bad_rg.wav'};
method = 'peri'; % peri, welch, fft, fftband
close all

for i = 1:length(filenames)
    file = filenames{i};
    [wave,Fs]=audioread(file);
    x=(wave(:,1)+wave(:,2))/2; %Downmix
    
    %     %Filter
    %     N   = 200;        % FIR filter order
    %     Fp  = 15000;       % 20 kHz passband-edge frequency
    %     Rp  = 0.00057565; % Corresponds to 0.01 dB peak-to-peak ripple
    %     Rst = 1e-4;       % Corresponds to 80 dB stopband attenuation
    %     eqnum = firceqrip(N,Fp/(Fs/2),[Rp Rst],'passedge'); % eqnum = vec of coeffs
    %     lowpassFIR = dsp.FIRFilter('Numerator',eqnum);
    %     x=lowpassFIR(x);
    
    x=normalizerms(x, -20);
    
    switch method
        case 'peri' % periodogram
            
            figure
            periodogram(x,rectwin(length(x)),length(x),Fs)
            
            xlim([0 20])
            ylim([-180 -20])
            
        case 'welch' % Welch's power spectral density estimate
            
            [pxx,f] = pwelch(x,2000,[],[],Fs); % Disclamer: I have zero idea about those parameters
            figure
            plot(f,10*log10(pxx)) % You can remove dB conversion
            xlabel('Frequency (Hz)')
            ylabel('Magnitude (dB)')
            
            xlim([0 20000])
            ylim([-105 -40])
            
        case 'fft' % FFT (amplitude spectrum)
            
            L=length(x);
            f = Fs*(0:(L/2))/L;
            Y = fft(x);
            P2 = abs(Y/L); % You can choose to not divide it by L, no difference
            P1 = P2(1:L/2+1);
            P1(2:end-1) = 2*P1(2:end-1);
            
            figure
            plot(f,P1) % you can convert it to dB too, make sure to use 'voltage'
            xlabel('f (Hz)')
            ylabel('|P1(f)|')
            
            ylim([0 0.004])
            xlim([0 20000])
            
        case 'fftband'
            
            L=length(x);
            f = Fs*(0:(L/2))/L;
            Y = fft(x); 
            P2 = abs(Y); % I don't divide by L here, because I'm going to use dB later
            P1 = P2(1:L/2+1);
            P1(2:end-1) = 2*P1(2:end-1);
            
            bars = 20;
%             % Linear spacing            
%             fvalues = linspace(0,20000,bars+1);            
%             fvalues = fvalues(2:end);
%             % Log10 spacing
%             fvalues = logspace(-1,log10(20000),bars);             
            % Copied from FB            
            fvalues=[50,69,94,129,176,241,331,453,620,850,1200,1600,2200,3000,4100,5600,7700,11000,14000,20000];     
            
            fband = zeros(length(fvalues),1);
            
            % Again, I don'tknow shoul you just add them together.
            for k = 1:length(f)
                for j = 1:length(fvalues)
                    if f(k) < fvalues(j)
                        fband(j) = fband(j)+P1(k);
                        break
                    end
                end
            end
            
            figure
            histogram('BinEdges',0:length(fvalues),'BinCounts',round(db(fband,'voltage'))')
            ylim([0 180])
            
            labels = [0 fvalues(2:2:length(fvalues))];
            xticklabels(num2cell(labels))            
            xlabel('f (Hz)')
            ylabel('dB')            
    end
    
    title(file)
    
end