Analysis of a musical sample after spectral filtering

We use Fourier and Butterworth filtering to break a short segment of Handel's Messiah into high, low and middle frequencies, which we then play back.

Contents

Load data and Fourier filter

  load handel % This music dataset in vector y is built into Matlab.
              % Note sampling frequency is Fs = 8192 Hz

  clf
  Nw = 2^16;  % Segment length for DFT (= 8 sec)- power of 2 for efficiency
  yw = y(1:Nw);  % Truncate signal to this length

  Mw = [0:(Nw/2 - 1) -Nw/2:-1];  % Harmonics
  fw = Fs*Mw/Nw; % "Frequencies'(inverse periods in Hz) of these harmonics

  % Find DFT indices corresponding to frequency ranges to zero out
  flt440 = (abs(fw) < 440)';
  flt880 = (abs(fw) < 880)';

  % DFT and high-pass filtering

  ywhat = fft(yw);

  yhi440hat = ywhat;
  yhi440hat(flt440,:) = 0;  % Zero out low frequencies < 440 Hz
  yhi440 = ifft(yhi440hat);

  yhi880hat = ywhat;
  yhi880hat(flt880,:) = 0;
  yhi880 = ifft(yhi880hat);

  ylo440 = yw - yhi440;
  yband  = yhi440 - yhi880;

  % Power in different spectral bands adds up to total power as it should
  power = mean(yw.^2)
  powerlo440 = mean(ylo440.^2)
  powerband = mean(yband.^2)
  powerhi880 = mean(yhi880.^2)

  % Plot short section of time series

  subplot(2,1,1)
  isamp = (8150):(8250);  % time index range to plot
  plot(isamp,y(isamp),'k',...
       isamp,ylo440(isamp),'r',...
       isamp,ylo440(isamp)+yband(isamp),'b')
  ylim([1.5*min(y(isamp)) 1.5*max(y(isamp))]);
  legend('y','<440 Hz','<880 Hz')
  title('Fourier filters')

power =

    0.0381


powerlo440 =

    0.0047


powerband =

    0.0188


powerhi880 =

    0.0146

Play the original time series

  p8 = audioplayer(y,Fs);
  playblocking(p8);

Low-pass filtered

  p8 = audioplayer(ylo440,Fs);
  playblocking(p8);

Band-pass filtered

  p8 = audioplayer(yband,Fs);
  playblocking(p8);

High-pass filtered

  p8 = audioplayer(yhi880,Fs);
  playblocking(p8);

Use 8th order Butterworth filter

  % Set up filters

  N = 8; % Filter order.
  W440 = 440/(Fs/2); % 440 Hz cutoff frequency (units of Nyquist freq Fs/2)
  W880 = 880/(Fs/2); % 880 Hz cutoff frequency (units of Nyquist freq Fs/2)
  [b440  a440]  = butter(N,W440); % Coefficients for 2N-order low-pass Butterworth filter
  [bband aband] = butter(N,[W440 W880]);
  [b880  a880]  = butter(N,W880,'high');

  % Filter the time series

  vlo440 = filter(b440,a440,y);
  vband = filter(bband,aband,y);
  vhi880 = filter(b880,a880,y);

  % Forward-backward filter the time series (removes phase delay)

  v2lo440 = filtfilt(b440,a440,y);
  v2band = filtfilt(bband,aband,y);
  v2hi880 = filtfilt(b880,a880,y);

  % Power of Butterworth-filtered signal in the three spectral bands
  % - compares well with the Fourier filters

  Bpowerlo440 = mean(vlo440.^2)
  Bpowerband = mean(vband.^2)
  Bpowerhi880 = mean(vhi880.^2)

  % Plot short section of filtered time series
  % -note how filtfilt removes the phase delay and makes filtered signal
  %  very similar to Fourier filter (but slightly attenuates some freqs.)

  subplot(2,1,2)
  isamp = (8150):(8250);
  plot(isamp,y(isamp),'k',...
       isamp,vlo440(isamp),'r--',...
       isamp,vlo440(isamp)+vband(isamp),'b--',...
       isamp,v2lo440(isamp),'r',...
       isamp,v2lo440(isamp)+v2band(isamp),'b')
  ylim([1.5*min(y(isamp)) 1.5*max(y(isamp))]);
  legend('y','<440 Hz','<880 Hz')
  title([int2str(N) 'th-order Butterworth filters'])
  text(0.96*isamp(1)+0.04*isamp(end),1.2*max(y(isamp)),'dash: 1-pass (filt)')
  text(0.96*isamp(1)+0.04*isamp(end),0.9*max(y(isamp)),'solid: 2-pass (filtfilt)')

Bpowerlo440 =

    0.0047


Bpowerband =

    0.0195


Bpowerhi880 =

    0.0146

Play low-pass Fourier filtered

  p8 = audioplayer(ylo440,Fs);
  playblocking(p8);

Play low-pass Butterworth filtered

  p8 = audioplayer(vlo440,Fs);
  playblocking(p8);

Play forward-backward low-pass Butterworth filtered

  p8 = audioplayer(v2lo440,Fs);
  playblocking(p8);

Play band-pass Fourier filtered

  p8 = audioplayer(yband,Fs);
  playblocking(p8);

Play band-pass Butterworth filtered

  p8 = audioplayer(vband,Fs);
  playblocking(p8);

Play high-pass Fourier filtered

  p8 = audioplayer(vhi880,Fs);
  playblocking(p8);

Play high-pass Butterworth filtered

  p8 = audioplayer(vhi880,Fs);
  playblocking(p8);

Published with MATLAB® R2012b