SPACE STUDY - DIRECTIONALITY

Contents

 Created by Juan-Pablo Caceres 30-Aug-2006.
 Last Revision: Under Subversion Control System
clear all
close all

Motivation

Filter design for impluse responses in order to remove the wall filtering effect, but keep the source position outside the office.

Impluse responses load and analysis

%-------------------------------------
load '~/yamaha/scripts/rooms/tokyo/tokyoIR'; %% Load impulse response (use imp_resp.m to generate)

Nspots = size(IRS,1);
Nmics = size(IRS,2);
ntaps = fs;

time = (0:ntaps-1)/fs*1000; %time vector (ms)

% Filter desing for high pass the impulses responses. This is used
% just to have a better fingerprint of where the direct path and
% the reflections are located in the impulses responses.
[b,a] = butter(6,500/44100*2,'high');

% Impulse response Adjacent room, Sound Filed Microphones
ir_adj = [IRS(2,5).ir IRS(2,6).ir IRS(2,7).ir IRS(2,8).ir];
ir_adjF = filter(b,a,[IRS(2,5).ir IRS(2,6).ir IRS(2,7).ir IRS(2,8).ir]);

maximp = max(abs(ir_adj));
plot(time,ir_adjF + ones(ntaps,4)*[0:3]'*[0 1 2 3]*.00003)
grid
xlabel('time (ms)')
legend({'dir01','dir02','dir03','dir04'})
title('Filtered Impulses Response (Tokyo Office), Adjacent Room')
datacursormode on %% To get points in the plot with the cursor
%ginput(2) %% Get cursor data to the command window

Spectral Smoothings

%---------------------
beta = 5.0;
nbins = 2^9;
ndirect = 5*nbins; %% Direct path length
[Ypsds1,f,Ypsd1] = specsmooth(ir_adj(1:ndirect,1),fs,nbins,beta);
[Ypsds2,f,Ypsd2] = specsmooth(ir_adj(1:ndirect,2),fs,nbins,beta);
[Ypsds3,f,Ypsd3] = specsmooth(ir_adj(1:ndirect,3),fs,nbins,beta);
[Ypsds4,f,Ypsd4] = specsmooth(ir_adj(1:ndirect,4),fs,nbins,beta);

% Plots
%-------
figure
semilogx(f/1000,10*log10(Ypsd1),'b')
hold on
semilogx(f/1000,10*log10(Ypsd2),'g')
semilogx(f/1000,10*log10(Ypsd3),'r')
semilogx(f/1000,10*log10(Ypsd4),'c')

semilogx(f/1000,10*log10(Ypsds1),'b','LineWidth',4)
semilogx(f/1000,10*log10(Ypsds2),'g','LineWidth',4)
semilogx(f/1000,10*log10(Ypsds3),'r','LineWidth',4)
semilogx(f/1000,10*log10(Ypsds4),'c','LineWidth',4)
legend({'PDS dir1';'PDS dir2';'PDS dir3';'PDS dir4';...
       'PDS smooth dir1';'PDS smooth dir2';'PDS smooth dir3';'PDS smooth dir4'})
ylabel('Power/Frequency (dB/Hz)')
xlabel('Frequency (kHz)')
title('PSD Impulse Response adjacent the room')
grid on
set(gca,'XLim',[0.05 22])

Filter pre-design (tweak)

-------------------------- Tweak the smooth PSD in the high frequencies in order to get a better minimum phase filter to invert We use PSD 3, since is the mic that is phasing the wall where the intruding sound is comming.

maxgain = 30; %max gain to apply in the tweaking =======================
Ypsds3N = Ypsds3/max(abs(Ypsds3));
Ypsds3NT = Ypsds3N + 10^(-maxgain/10); %Normalize and tweak
Ypsd3N = Ypsd3/max(abs(Ypsd3));

Filter design (minimum phase)

%------------------------------
ys = real(ifft(exp(conj(hilbert(1/2*log([Ypsds3NT; flipud(Ypsds3NT(2:nbins))]))))));
nu = .3; %KHz frequency in the middle of the band (magnifying glass)
order = 8;
[bw, aw] = pronyw(ys(1:nbins), order, [nu*2000/fs 0.5]);
tfhatw = freqz(bw, aw, pi*[0:nbins]'/nbins);
temp = fft(ys(1:nbins), 2*nbins);


% Plots
%-------
figure
semilogx(f/1000,10*log10(Ypsd3N),'r')
hold on
semilogx(f/1000,10*log10(Ypsds3N),'r','LineWidth',4)
semilogx(f/1000,10*log10(Ypsds3NT),'k--','LineWidth',2)
semilogx(f/1000, 20*log10(abs(tfhatw)),'g','LineWidth',2);
legend({'PDS dir3';'PDS smooth dir3';'PDS smooth Tweak dir3';'IIR Filter'})
ylabel('Power/Frequency (dB/Hz)')
xlabel('Frequency (kHz)')
title('PSD Impulse Response adjacent the room, mic dir03')
grid on
set(gca,'XLim',[0.05 22])

Inverse filter Design

%-----------------------
binv = aw*sum(bw)/sum(aw); %% normalized
ainv = bw;
% ir_adjIF is equalized impulse responses to be used to remove the wall effect.
ir_adjIF = filter(binv,ainv,ir_adj);


figure
subplot(211)
specgram(ir_adj(:),256,fs/1000)
grid
colorbar
set(gca,'YLim',[0 5])
subplot(212)
specgram(ir_adjIF(:),256,fs/1000)
grid
set(gca,'YLim',[0 5])
colorbar

Sound

%-------
% Filter to remove the high (noise) frequencies in the impulse response
[b,a] = butter(2,7000/44100*2);
%freqz(b,a)
wavwrite(0.9*ir_adj(:)/max(abs(ir_adj(:))),fs,16,'original.wav')
wavwrite(0.9*filter(b,a,ir_adjIF(:)/max(abs(ir_adjIF(:)))),fs,16,'eq.wav')
wavwrite(0.9*ir_adjIF(:)/max(abs(ir_adjIF(:))),fs,16,'eqnf.wav')

ir_2 = [IRS(1,5).ir IRS(1,6).ir IRS(1,7).ir IRS(1,8).ir];
wavwrite(0.9*ir_2(:)/max(abs(ir_2(:))),fs,16,'inside.wav')

Save IR data into a struct

%----------------------------

for i = 1:4
  IRSnowall(1,i).ir = filter(b,a,ir_adjIF(:,i));
  IRSnowall(1,i).spot = IRS(2,i+4).spot;
  IRSnowall(1,i).mic = IRS(2,i+4).mic;
end

Testing

%---------

Generation of the convolution of the conversations with the measures impulse responses.

%Select output path where the wav files will be saved.
in_path = '~/yamaha/recordings/experiments/conversations/';
in_filenames = {'2people_e.wav'};
out_path = '~/yamaha/recordings/space_study/CONVconversations/';
load tokyoIR;
ntaps = length(IRS(1,1).ir)/2;
%convgen(IRS,ntaps,fs,in_path,out_path,in_filenames) %Convolution


out_path = '~/yamaha/recordings/space_study/CONVconversationsNW/';
%convgen(IRSnowall,ntaps,fs,in_path,out_path,in_filenames) %Convolution
Error using ==> load
Unable to read file tokyoIR: No such file or directory.