0001 function [c,qq,ff,tt,po]=v_modspect(s,fs,m,nf,nq,p)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092
0093
0094
0095
0096
0097
0098
0099
0100
0101
0102
0103
0104
0105
0106
0107
0108
0109
0110
0111
0112
0113
0114
0115
0116
0117
0118
0119
0120
0121
0122
0123
0124
0125
0126
0127
0128
0129
0130
0131
0132
0133 persistent PP
0134 if isempty(PP)
0135
0136
0137
0138
0139 PP.logr=120;
0140 PP.tagc=0;
0141 PP.tagt=2;
0142 PP.tinc=0.25e-3;
0143 PP.tovl=16;
0144 PP.tpad=30e-3;
0145
0146 PP.trnd=1;
0147 PP.twin='m';
0148 PP.fpow='p';
0149 PP.fmin=40;
0150 PP.fmax=10000;
0151 PP.fwin='t';
0152 PP.fbin=0.1;
0153 PP.ppow='a';
0154 PP.pinc=16e-3;
0155 PP.povl=2;
0156 PP.pwin='m';
0157 PP.ppad=200e-3;
0158
0159 PP.prnd=1;
0160 PP.mpow='p';
0161 PP.mwin='t';
0162 PP.mbin=15;
0163 PP.mmin=50;
0164 PP.mmax=400;
0165 PP.qpow='a';
0166 PP.qdcq=0;
0167 PP.qdcp=0;
0168 PP.dzcq=0;
0169 PP.dzcp=0;
0170 PP.dncp=30;
0171 PP.dncq=15;
0172 PP.ddet=0;
0173 PP.ddnt=2;
0174 PP.ddep=0;
0175 PP.ddnp=1;
0176 PP.unit=0;
0177 PP.dplt=0;
0178
0179 PP.cent=0.2;
0180 PP.cchk=0.2;
0181 end
0182 po=PP;
0183 switch nargin
0184 case 0
0185
0186 nn=sort(fieldnames(PP));
0187 cnn=char(nn);
0188 fprintf('%d Voicebox parameters:\n',length(nn));
0189
0190 for i=1:length(nn);
0191 if ischar(PP.(nn{i}))
0192 fmt=' %s = %s\n';
0193 else
0194 fmt=' %s = %g\n';
0195 end
0196 fprintf(fmt,cnn(i,:),PP.(nn{i}));
0197 end
0198 return
0199 case 1
0200 error('no sample frequency specified');
0201 case 2
0202 p=[];
0203 case 3
0204 if isstruct(m)
0205 p=m;
0206 m='';
0207 else
0208 p=[];
0209 end
0210 nf=[];
0211 nq=[];
0212 case 4
0213 if isstruct(nf)
0214 p=nf;
0215 nf=[];
0216 else
0217 p=[];
0218 end
0219 nq=[];
0220 case 5
0221 if isstruct(nq)
0222 p=nq;
0223 nq=[];
0224 else
0225 p=[];
0226 end
0227 end
0228 if isstruct(p)
0229 nn=fieldnames(p);
0230 for i=1:length(nn)
0231 if isfield(po,nn{i})
0232 po.(nn{i})=p.(nn{i});
0233 end
0234 end
0235 end
0236
0237
0238
0239 for i=1:length(m)
0240 switch m(i)
0241 case 'g'
0242 po.tagc=1;
0243 case {'r','n','m'}
0244 po.twin=m(i);
0245 po.pwin=m(i);
0246 case {'T','N','M'}
0247 po.fwin=lower(m(i));
0248 po.mwin=po.fwin;
0249 case {'a','c'}
0250 po.fpow=m(i);
0251 case {'p','l'}
0252 po.ppow=m(i);
0253 case 'A'
0254 po.mpow=lower(m(i));
0255 case {'P','L'}
0256 po.qpow=lower(m(i));
0257 case 'f'
0258 po.qdcq=1;
0259 case 'F'
0260 po.qdcp=1;
0261 case 'z'
0262 po.dzcq=1;
0263 case 'Z'
0264 po.dzcp=1;
0265 case 'd'
0266 po.ddet=1;
0267 case 'D'
0268 po.ddep=1;
0269 case 'u'
0270 po.unit=1;
0271 case 's'
0272 po.dplt=bitor(po.dplt,31+512);
0273 case 'S'
0274 po.dplt=bitor(po.dplt,1023);
0275 end
0276 end
0277
0278 if ~nargout
0279 po.dplt=bitor(po.dplt,4);
0280 end
0281 if ~isempty(nf)
0282 po.dncp=nf;
0283 end
0284 if ~isempty(nq)
0285 po.dncq=nq;
0286 end
0287
0288 lnf=length(nf);
0289 if lnf>=1
0290 po.fbin=nf(1);
0291 if lnf>=2
0292 po.fmin=nf(2);
0293 if lnf>=3
0294 po.fmax=nf(3);
0295 if lnf>=4
0296 po.dncp=nf(4);
0297 end
0298 end
0299 end
0300 end
0301 lnq=length(nq);
0302 if lnq>=1
0303 po.mbin=nq(1);
0304 if lnq>=2
0305 po.mmin=nq(2);
0306 if lnq>=3
0307 po.mmax=nq(3);
0308 if lnq>=4
0309 po.dncq=nq(4);
0310 end
0311 end
0312 end
0313 end
0314
0315
0316
0317 po.ddep=po.ddep & (po.qdcp==0);
0318 po.dzcq=po.dzcq | (po.qdcq==0);
0319 po.dzcp=po.dzcp | (po.qdcp==0);
0320
0321 logth=10^(-po.logr/10);
0322 ns=length(s);
0323 axhandle=[];
0324
0325
0326
0327 if po.tagc
0328 tau=po.tagt*fs;
0329 ax2=[1 -exp(-1/tau)];
0330 bx2=sum(ax2);
0331 sm=mean(s(1:min(ns,round(tau))).^2);
0332 if sm>0
0333 s=s./sqrt(filter(bx2,ax2,s.^2,-ax2(2)*sm));
0334 end
0335 end
0336 if bitand(po.dplt,32)
0337 if bitand(po.dplt,1)
0338 figure;
0339 end
0340 plot((1:ns)/fs,s);
0341 axhandle(end+1)=gca;
0342 title('Input after AGC');
0343 xlabel('Time (s)');
0344 end
0345
0346
0347
0348 ni=ceil(po.tinc*fs);
0349 nw=ni*po.tovl;
0350 nf=ceil(max(nw,po.tpad*fs));
0351 if po.trnd
0352 nf=2^nextpow2(nf);
0353 else
0354 nf=round(nf);
0355 end
0356 switch po.twin
0357 case 'm'
0358
0359 w=0.54-0.46*cos(2*pi*(0:nw-1)/nw);
0360 case 'n'
0361
0362 w=0.5-0.5*cos(2*pi*(1:nw-1)/nw);
0363 case 'r'
0364 w=ones(nw,1);
0365 end
0366 fx=v_rfft(v_enframe(s,w,ni),nf,2);
0367 [nft,nff]=size(fx);
0368 if bitand(po.dplt,64)
0369 if bitand(po.dplt,1)
0370 figure;
0371 end
0372 fp=fx.*conj(fx);
0373 imagesc(((0:nft-1)*ni+(nw+1)/2)/fs,(0:nff-1)*fs*0.001/nf,10*log10(max(fp,max(fp(:))*1e-6))');
0374 axhandle(end+1)=gca;
0375 hanc= colorbar;
0376 set(get(hanc,'ylabel'),'String','dB');
0377 axis('xy');
0378 title('Input Spectrogram');
0379 xlabel('Time (s)');
0380 ylabel('Frequency (kHz)');
0381 end
0382 [mbk,ffm]=v_melbankm(po.fbin,nf,fs,po.fmin/fs,min(po.fmax/fs+(po.fmax==0),0.5),po.fwin);
0383 switch [po.fpow po.ppow]
0384 case 'aa'
0385 px=abs(fx)*mbk';
0386 case {'ap','al'}
0387 px=(abs(fx)*mbk').^2;
0388 case 'pa'
0389 px=sqrt((fx.*conj(fx))*mbk');
0390 case {'pp','pl'}
0391 px=(fx.*conj(fx))*mbk';
0392 case 'ca'
0393 px=abs(fx*mbk');
0394 case {'cp','cl'}
0395 px=fx*mbk';
0396 px=px.*conj(px);
0397 end
0398 if po.ppow=='l'
0399 px=log(max((px),max(px(:))*logth));
0400 end
0401 if bitand(po.dplt,128)
0402 if bitand(po.dplt,1)
0403 figure;
0404 end
0405 switch po.ppow
0406 case 'a'
0407 imagesc(((0:nft-1)*ni+(nw+1)/2)/fs,ffm/1000,20*log10(max(px,max(px(:))*1e-3))');
0408 case 'p'
0409 imagesc(((0:nft-1)*ni+(nw+1)/2)/fs,ffm/1000,10*log10(max(px,max(px(:))*1e-6))');
0410 case 'l'
0411 imagesc(((0:nft-1)*ni+(nw+1)/2)/fs,ffm/1000,max(px*10/log(10),max(px(:))*10/log(10)-60)');
0412 end
0413 axhandle(end+1)=gca;
0414 hanc= colorbar;
0415 set(get(hanc,'ylabel'),'String','dB');
0416 axis('xy');
0417 title('Mel Spectrogram');
0418 xlabel('Time (s)');
0419 ylabel('Frequency (k-mel)');
0420 end
0421 npf=length(ffm);
0422 mni=ceil(po.pinc*fs/ni);
0423 mnw=mni*po.povl;
0424 mnf=ceil(max(mnw,po.ppad*fs/ni));
0425 if po.prnd
0426 mnf=2^nextpow2(mnf);
0427 else
0428 mnf=round(mnf);
0429 end
0430 nmt=fix((nft-mnw+mni)/mni);
0431 ix=repmat((1:mnw)',1,nmt)+repmat((0:nmt-1)*mni,mnw,1);
0432 mx=v_rfft(reshape(px(ix(:),:),mnw,nmt*npf),mnf,1);
0433
0434 [qbk,qqm]=v_melbankm(po.mbin,mnf,100*fs/ni,po.mmin*ni/fs,min(po.mmax*ni/fs+(po.mmax==0),0.5),po.mwin);
0435 nqq=length(qqm);
0436 switch po.mpow
0437 case 'a'
0438 qx=reshape(qbk*abs(mx),nqq,nmt,npf);
0439 case 'p'
0440 qx=reshape(qbk*(mx.*conj(mx)),nqq,nmt,npf);
0441 end
0442 switch [po.mpow po.qpow]
0443 case 'aa'
0444 qx=reshape(qbk*abs(mx),nqq,nmt,npf);
0445 case {'ap','al'}
0446 qx=reshape(qbk*abs(mx),nqq,nmt,npf).^2;
0447 case 'pa'
0448 qx=sqrt(reshape(qbk*(mx.*conj(mx)),nqq,nmt,npf));
0449 case {'pp','pl'}
0450 qx=reshape(qbk*(mx.*conj(mx)),nqq,nmt,npf);
0451 end
0452 if po.qpow=='l'
0453 qx=log(max((qx),max(qx(:))*logth));
0454 end
0455 tt=((1+mnw*ni+nw-ni)/2+(0:nmt-1)*mni*ni)/fs;
0456
0457 if bitand(po.dplt,256) && (~bitand(po.dplt,4) || po.qdcq>0 || po.qdcp>0)
0458 if bitand(po.dplt,1)
0459 figure;
0460 end
0461 switch po.qpow
0462 case 'a'
0463 dqx=20*log10(max(qx,max(qx(:))*1e-3));
0464 case 'p'
0465 dqx=10*log10(max(qx,max(qx(:))*1e-6));
0466 case 'l'
0467 dqx=max(qx*10/log(10),max(qx(:))*10/log(10)-60);
0468 end
0469 dqx(end+1,:,:)=max(dqx(:));
0470 dqx=reshape(permute(dqx,[2,1,3]),nmt,(nqq+1)*npf);
0471 ffq=ffm(1)+((0:(nqq+1)*npf)-(nqq-1)/2)*(ffm(2)-ffm(1))/(nqq+1);
0472 imagesc(tt,ffq/1000,dqx');
0473 axhandle(end+1)=gca;
0474 hanc= colorbar;
0475 set(get(hanc,'ylabel'),'String','dB');
0476 axis('xy');
0477 title('Modulation Spectrogram');
0478 xlabel('Time (s)');
0479 ylabel('Frequency (k-mel)');
0480 end
0481 ndq=nqq;
0482 if po.qdcq
0483 dx=reshape(v_rdct(reshape(qx,nqq,nmt*npf)),nqq,nmt,npf);
0484 ndq=min(ndq,po.dncq+1);
0485 else
0486 dx=qx;
0487 end
0488 ndf=npf;
0489 if po.qdcp
0490 dx=permute(reshape(v_rdct(reshape(permute(qx,[3,1,2]),npf,nqq*nmt)),npf,nqq,nmt),[2 3 1]);
0491 ndf=min(ndf,po.dncp+1);
0492 elseif po.ddep
0493 nv=po.ddnp;
0494 vv=(-nv:nv)*-3/(nv*(nv+1)*(2*nv+1)*(ffm(2)-ffm(1)));
0495 dxdp=filter(vv,1,dx,[],3);
0496 dxdp=dxdp(:,:,[repmat(2*nv+1,1,nv) 2*nv+1:npf repmat(npf,1,nv)]);
0497 end
0498 if po.ddet
0499 nv=po.ddnt;
0500 vv=(-nv:nv)*-3/(nv*(nv+1)*(2*nv+1)*(tt(2)-tt(1)));
0501 dxdt=filter(vv,1,dx,[],2);
0502 dxdt=dxdt(:,[repmat(2*nv+1,1,nv) 2*nv+1:nmt repmat(nmt,1,nv)],:);
0503 end
0504 nqj=ndq-(po.dzcq==0);
0505 nqk=nqj+ndq*(po.ddep+po.ddet);
0506 npk=ndf-(po.dzcp==0);
0507 c=zeros(nqk,npk,nmt);
0508 c(1:nqj,:,:)=permute(dx(1+ndq-nqj:ndq,:,1+ndf-npk:ndf),[1,3,2]);
0509 if bitand(po.dplt,512)
0510 if bitand(po.dplt,1)
0511 figure;
0512 end
0513 ffx=repmat(v_mel2frq(ffm(:)),1,nqq)/1000;
0514 qqx=repmat(mel2frq(qqm(:)')/100,npf,1);
0515 plot(qqx,ffx,'xb');
0516 axis([[qqx(1) qqx(end)]*[1.1 -0.1; -0.1 1.1] [ffx(1) ffx(end)]*[1.1 -0.1; -0.1 1.1]]);
0517 title('Frequency bin centres');
0518 xlabel(sprintf('%.1f : %.1f : %.1f mod-mel = Modulation Frequency (Hz)',qqm(1)/100,(qqm(2)-qqm(1))/100,qqm(end)/100));
0519 ylabel(sprintf('%.0f : %.0f : %.0f mel = Frequency (kHz)',ffm(1),ffm(2)-ffm(1),ffm(end)));
0520 end
0521 if bitand(po.dplt,4)
0522 if bitand(po.dplt,1)
0523 figure;
0524 end
0525 dqx=dx(1+ndq-nqj:ndq,:,1+ndf-npk:ndf);
0526 dqxmx=max(abs(dqx(:)));
0527 dqxge=dqx>=0;
0528 dbfact=2-(po.qpow=='p');
0529 dqxa=max(abs(dqx),dqxmx*10^(-6/dbfact));
0530 dqxmn=min(dqxa(:));
0531 dblab='';
0532 if(mean(dqxa(:)>dqxmn+po.cent*(dqxmx-dqxmn)))<po.cchk
0533 dboff=abs(round(dbfact*10*log10(dqxmn)));
0534 dblab='{\pm}dB';
0535 if ~all(dqxge(:))
0536 dqxa=dqxa/dqxmn;
0537 if dqxmn>1 && dboff~=0
0538 dblab=sprintf('{\\pm}dB - %d',dboff);
0539 else
0540 dblab=sprintf('{\\pm}dB + %d',dboff);
0541 end
0542 else
0543 dblab='dB';
0544 end
0545 dqx=dbfact*10*log10(dqxa).*(2*dqxge-1);
0546 end
0547 dqx(end+1,:,:)=max(dqx(:));
0548 dqx=reshape(permute(dqx,[2,1,3]),nmt,(nqj+1)*npk);
0549 ffq=ffm(1)+((0:(nqj+1)*npf)-(nqj-1)/2)*(ffm(2)-ffm(1))/(nqj+1);
0550 imagesc(tt,ffq/1000,dqx');
0551 axhandle(end+1)=gca;
0552 hanc= colorbar;
0553 set(get(hanc,'ylabel'),'String',dblab);
0554 axis('xy');
0555 if po.qdcq
0556 title('Modulation spectrum DCT');
0557 else
0558 title('Modulation spectrum');
0559 end
0560 xlabel('Time (s)');
0561 if po.qdcp
0562 ylabel('Quefrency (k-mel^{-1})');
0563 else
0564 ylabel('Frequency (k-mel)');
0565 end
0566 end
0567
0568 if po.ddep
0569 c(nqj+1:nqj+ndq,:,:)=permute(dxdp(1:ndq,:,1+ndf-npk:ndf),[1,3,2]);
0570 if bitand(po.dplt,8)
0571 if bitand(po.dplt,1)
0572 figure;
0573 end
0574 dqx=dxdp(1:ndq,:,1+ndf-npk:ndf);
0575 dqxmx=max(abs(dqx(:)));
0576 dqxge=dqx>=0;
0577 dbfact=2-(po.qpow=='p');
0578 dqxa=max(abs(dqx),dqxmx*10^(-6/dbfact));
0579 dqxmn=min(dqxa(:));
0580 dblab='';
0581 if(mean(dqxa(:)>dqxmn+po.cent*(dqxmx-dqxmn)))<po.cchk
0582 dboff=abs(round(dbfact*10*log10(dqxmn)));
0583 dblab='{\pm}dB';
0584 if ~all(dqxge(:))
0585 dqxa=dqxa/dqxmn;
0586 if dqxmn>1 && dboff~=0
0587 dblab=sprintf('{\\pm}dB - %d',dboff);
0588 else
0589 dblab=sprintf('{\\pm}dB + %d',dboff);
0590 end
0591 else
0592 dblab='dB';
0593 end
0594 dqx=dbfact*10*log10(dqxa).*(2*dqxge-1);
0595 end
0596 dqx(end+1,:,:)=max(dqx(:));
0597 dqx=reshape(permute(dqx,[2,1,3]),nmt,(ndq+1)*npk);
0598 ffq=ffm(1)+((0:(ndq+1)*npf)-(ndq-1)/2)*(ffm(2)-ffm(1))/(ndq+1);
0599 imagesc(tt,ffq/1000,dqx');
0600 axhandle(end+1)=gca;
0601 hanc= colorbar;
0602 set(get(hanc,'ylabel'),'String',dblab);
0603 axis('xy');
0604 if po.qdcq
0605 title('Modulation spectrum DCT = freq derivative');
0606 else
0607 title('Modulation spectrum = freq derivative');
0608 end
0609 xlabel('Time (s)');
0610 if po.qdcp
0611 ylabel('Quefrency (k-mel^{-1})');
0612 else
0613 ylabel('Frequency (k-mel)');
0614 end
0615 end
0616 end
0617 if po.ddet
0618 c(nqk-ndq+1:nqk,:,:)=permute(dxdt(1:ndq,:,1+ndf-npk:ndf),[1,3,2]);
0619 if bitand(po.dplt,16)
0620 if bitand(po.dplt,1)
0621 figure;
0622 end
0623 dqx=dxdt(1:ndq,:,1+ndf-npk:ndf);
0624 dqxmx=max(abs(dqx(:)));
0625 dqxge=dqx>=0;
0626 dbfact=2-(po.qpow=='p');
0627 dqxa=max(abs(dqx),dqxmx*10^(-6/dbfact));
0628 dqxmn=min(dqxa(:));
0629 dblab='';
0630 if(mean(dqxa(:)>dqxmn+po.cent*(dqxmx-dqxmn)))<po.cchk
0631 dboff=abs(round(dbfact*10*log10(dqxmn)));
0632 dblab='{\pm}dB';
0633 if ~all(dqxge(:))
0634 dqxa=dqxa/dqxmn;
0635 if dqxmn>1 && dboff~=0
0636 dblab=sprintf('{\\pm}dB - %d',dboff);
0637 else
0638 dblab=sprintf('{\\pm}dB + %d',dboff);
0639 end
0640 else
0641 dblab='dB';
0642 end
0643 dqx=dbfact*10*log10(dqxa).*(2*dqxge-1);
0644 end
0645 dqx(end+1,:,:)=max(dqx(:));
0646 dqx=reshape(permute(dqx,[2,1,3]),nmt,(ndq+1)*npk);
0647 ffq=ffm(1)+((0:(nqj+1)*npf)-(nqj-1)/2)*(ffm(2)-ffm(1))/(nqj+1);
0648 imagesc(tt,ffq/1000,dqx');
0649 axhandle(end+1)=gca;
0650 hanc= colorbar;
0651 set(get(hanc,'ylabel'),'String',dblab);
0652 axis('xy');
0653 if po.qdcq
0654 title('Modulation spectrum DCT = time derivative');
0655 else
0656 title('Modulation spectrum = time derivative');
0657 end
0658 xlabel('Time (s)');
0659 if po.qdcp
0660 ylabel('Quefrency (k-mel^{-1})');
0661 else
0662 ylabel('Frequency (k-mel)');
0663 end
0664 end
0665 end
0666 if po.unit
0667 ff=ffm;
0668 qq=qqm/100;
0669 else
0670 ff=v_mel2frq(ffm);
0671 qq=v_mel2frq(qqm)/100;
0672 end
0673 if length(axhandle)>1
0674 if ~bitand(po.dplt,32+64)
0675 linkaxes(axhandle)
0676 else
0677 linkaxes(axhandle,'x')
0678 end
0679 end
0680
0681