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.