www.pudn.com > psychoacoustic.rar > psychoacoustic.m, change:2011-07-11,size:12763b


 
LENGTH=1024 
 
%%数据表格 
%%%%%%%%汉宁窗%%%%%	 
for	i=0:1023 
	window(i+1)=1.632993*0.5*(1-cos(2*pi*i/LENGTH)); 
end 
 
%%%%%%%%%抽选后的索引值对应FFT的k%%%%%%长度132,P307,K:0-511,INDEX:0-130,size(index_k)=131; 
index_k=[0,1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32,33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48,50, 52, 54, 56, 58, 60, 62, 64, 66, 68, 70, 72, 74, 76, 78, 80, 82, 84, 86, 88, 90, 92, 94, 96,100,104,108,112,116,120,124,128,132,136,140,144,148,152,156,160,164,168,172,176,180,184,188,192,200,208,216,224,232,240,248,256,264,272,280,288,296,304,312,320,328,336,344,352,360,368,376,384,392,400,408,416,424,432,440,448,456,464]; 
%未被抽选到的序号处值为0,size(decimate)=465 
decimate=[0,1, 2, 3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29,30, 31, 32,33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48,0, 49, 0, 50, 0, 51, 0, 52, 0, 53, 0, 54, 0, 55, 0, 56,0, 57, 0, 58, 0, 59, 0, 60, 0, 61, 0, 62, 0, 63, 0, 64,0, 65, 0, 66, 0, 67, 0, 68, 0, 69,  0, 70, 0, 71, 0, 72,0, 0, 0, 73, 0, 0, 0, 74, 0, 0, 0, 75, 0, 0, 0, 76,0, 0, 0, 77, 0, 0, 0, 78, 0,0, 0, 79, 0, 0, 0, 80,0, 0, 0, 81, 0, 0, 0, 82, 0, 0, 0, 83, 0, 0, 0, 84,0, 0, 0, 85, 0, 0, 0,  86, 0, 0, 0, 87, 0, 0, 0, 88,0, 0, 0, 89, 0, 0, 0, 90, 0, 0, 0, 91, 0, 0, 0, 92,0, 0, 0, 93, 0, 0, 0, 94, 0, 0, 0, 95, 0, 0, 0, 96,0, 0, 0,0, 0, 0, 0, 97, 0, 0,0,0,0,0,0,98,0,0,0,0,0,0,0, 99,0,0,0,0,0,0,0,100,0,0,0,0,0,0,0,101,0,0,0,0,0,0,0,102,0,0,0,0,0,0,0,103,0,0,0,0,0,0,0,104,0,0,0,0,0,0,0,105,0,0,0,0,0,0,0,106,0,0,0,0,0,0,0,107,0,0,0,0,0,0,0,108,0,0,0,0,0,0,0,109,0,0,0,0,0,0,0,110,0,0,0,0,0,0,0,111,0,0,0,0,0,0,0,112,0,0,0,0,0,0,0,113,0,0,0,0,0,0,0,114,0,0,0,0,0,0,0,115,0,0,0,0,0,0,0,116,0,0,0,0,0,0,0,117,0,0,0,0,0,0,0,118,0,0,0,0,0,0,0,119,0,0,0,0,0,0,0,120,0,0,0,0,0,0,0,121,0,0,0,0,0,0,0,122,0,0,0,0,0,0,0,123,0,0,0,0,0,0,0,124,0,0,0,0,0,0,0,125,0,0,0,0,0,0,0,126,0,0,0,0,0,0,0,127,0,0,0,0,0,0,0,128,0,0,0,0,0,0,0,129,0,0,0,0,0,0,0,130];		 
		%长度465 
%FFT的k对应的抽选后的索引值 
for i=464:-1:0 
	if decimate(i+1)~=0 
		k_index(i+1)=decimate(i+1); 
	else 
		k_index(i+1)=k_index(i+1+1); 
	end 
end 
for i=465:511 
	k_index(i+1)=k_index(464+1); 
end 
 
%%%%%%%%%%27个临界带的BOTTOM,长度29,CRITBAND:0-26,INDEX:0-130,k:0-465,P317%%%%%%% 
boundindex=[0,1,2,3,4,6,8,11,14,17,20,23,27,31,36,41,47,52,57,63,70,75,80,86,93,100,106,118,131];%29 
boundk=[0,1,2,3,4,6,8,11,14,17,20,23,27,31,36,41,47,55,65,77,91,105,125,149,177,217,265,361,465];%29 
%51个半个临界带的BOTTOM, 长度52,p317,CRITBAND:0-50,K:0-464 
halfbarkk=[0,1,2,3,4,5,6,8,9,11,12,14,15,17,18,20,21,23,25,27,29,31,33,36,38,41,44,47,51,55,60,65,71,77,83,91,98,105,115,125,137,149,163,177,197,217,241,265,313,361,413,465]; 
 
%26个临界带的BARK值,长度131,P307,CRIBAND:0-26,INDEX:0-130 
CBRz=[0,0.425, 0.850, 1.273, 1.694, 2.112, 2.525, 2.934, 3.337, 3.733, 4.124,4.507, 4.882, 5.249, 5.608, 5.959, 6.301, 6.634, 6.959, 7.274, 7.581, 7.879, 8.169, 8.450, 8.723, 8.987, 9.244, 9.493, 9.734, 9.968,10.195,11.416,10.629,10.836,11.037,11.232,11.421,11.605,11.783,11.957,12.125,12.289,12.448,12.603,12.753,12.900,13.042,13.181,13.317,13.578,13.826,14.062,14.288,14.504,14.711,14.909,15.100,15.284,15.460,15.631,15.796,15.955,16.110,16.260,16.406,16.547,16.685,16.820,16.951,17.079,17.205,17.327,17.447,17.680,17.905,18.121,18.331,18.534,18.731,18.922,19.108,19.289,19.464,19.635,19.801,19.963,20.120,20.273,20.421,20.565,20.705,20.840,20.972,21.099,21.222,21.342,21.457,21.677,21.882,22.074,22.253,22.420,22.576,22.721,22.857,22.984,23.102,23.213,23.317,23.415,23.506,23.592,23.673,23.749,23.821,23.888,23.952,24.013,24.070,24.125,24.176,24.225,24.271,24.316,24.358,24.398,24.436,24.473,24.508,24.542,24.574];%131 
for i=1:464 
	if decimate(i+1)~=0 
		CBRzk(i+1)=CBRz(decimate(i+1)+1); 
	else 
		CBRzk(i+1)=CBRzk(i); 
	end 
end 
for i=465:511 
	CBRzk(i+1)=CBRzk(464+1); 
end 
 
%%被抽选到的索引处的静态掩蔽阈值,长度131,P307,CRIBAND:0-26,INDEX:0-130%%%%%%%%%%%%%%%% 
LTqindex=[0, 45.05, 25.87, 18.70, 14.85, 12.41, 10.72,  9.47,  8.05,  7.73,  7.10, 6.56,  6.11,  5.72,  5.37,  5.07,	4.79,  4.55,  4.32,  4.11,  3.92, 3.74,  3.57,  3.40,  3.25,  3.10,2.95,  2.81,  2.67,  2.53,  2.39, 2.25,  2.11,  1.97,  1.83,  1.68,	1.53,  1.38,  1.23,  1.07,  0.90, 0.74,  0.56,  0.39,  0.21,  0.02, -0.17, -0.36, -0.56, -0.96, -1.37, -1.79, -2.21, -2.63, -3.03, -3.41, -3.77, -4.09, -4.37, -4.60, -4.78, -4.91, -4.97, -4.98, -4.92, -4.81, -4.65, -4.43, -4.17, -3.87, -3.54, -3.19, -2.82, -2.06, -1.32, -0.64, -0.04,  0.47,  0.89,  1.23,  1.51,1.74,  1.93,  2.11,  2.28,  2.45,	2.63,  2.82,  3.03,  3.25,  3.49, 3.74,  4.02,  4.32,  4.64,  4.98,	5.35,  6.15,  7.07,  8.10,  9.25,10.54, 11.97, 13.56, 15.30, 17.23, 19.33, 21.64, 24.15, 26.88, 29.84, 33.04, 36.51, 40.24, 44.26, 48.58, 53.21, 58.17, 63.48, 68.00, 68.00, 68.00, 68.00, 68.00, 68.00, 68.00, 68.00, 68.00, 68.00, 68.00, 68.00]; 
%LTqindex=LTqindex+12;		 
%DFT的每个k处对应的静态掩蔽阈值 
for i=1:464 
	if decimate(i+1)~=0 
		LTqk(i+1)=LTqindex(decimate(i+1)+1); 
	else 
		LTqk(i+1)=LTqk(i); 
	end 
end 
for i=465:511 
	LTqk(i+1)=LTqk(464+1); 
end 
%%%%%%%%%%%%%%数据表格结束 
 
 
 
  faudin = fopen('E:\MATLAB\mywork\yang\mingyun.wav','rb'); 
	x=fread(faudin, 44,'char'); 
  fseek(faudin,100*1024,'bof'); 
%%%%%%%%%%%读一帧音频%%%%%%%%%%%%%%%% 
	x=fread(faudin,LENGTH,'int16') ;  %列矢量 
	figure; plot([0:1023],x);	xlabel('sample'); ylabel('amplitute');%figure1 
	 
	 
%%%%%%%%%%%%%%%%%应用心理声学模型%%%%%%%% 
	%fft的序号k:0-1023;	静态掩蔽域值和临界带率的序号index:1-130; 
	%临界带序号:critband:0:26 
	 
	%%%%%%%(1)fft()%%%%%%%%%% 
	Xabspow=(abs(fft(window.*x'))/LENGTH).^2; 
	X=10*log10(Xabspow);		m=max(X);		normal=96-m;	X=X+normal; 
					tmp=10.^(normal/10);	Xabspow=Xabspow*tmp; 
  figure;	plot([0:511]*44.1/1024,X(1:512),'r');	xlabel('frequency(kHz)'); ylabel('SPL(dB)');hold on;	%figuer2			 
	 
		%%%%%%%(2)TonalPressure(),用fft序号%%%%%%%%% 
	Xtm=zeros(1,512);	Xnm=Xtm;			tonal_flag=Xtm; 
  	for k=2:465-1 
		if  X(k+1)>X(k-1+1) & X(k+1)>=X(k+1+1) 
	    		tonal_flag(k+1)=2; 
			j=2; 
			if( k>=63 ) 
				j=3; 
			end 
			if( k>=127 ) 
				j=6; 
			end 
			if( k>=255 ) 
				j=12; 
			end 
			for ii=2:j 
		      		if X(k+1)-X(k-ii+1)<7 | X(k+1)-X(k+ii+1)<7 
			    		tonal_flag(k+1)=0; 
			    		ii=j+1; 
		      		end 
		      	end 
		end 
		if tonal_flag(k+1)==2  
      			Xtm(k+1)=10*log10( Xabspow(k-1+1)+Xabspow(k+1)+Xabspow(k+1+1) );	%DB 
		      	for ii=0:j 
		      		Xabspow(k-ii+1)=0; 
		      		Xabspow(k+ii+1)=0; 
		      	end 
	        	k=k+j+1; 
    end 
  end 
figure; h1=stem([0:511]*44.1/1024,Xtm);		%stem([0:511],tonal_flag*10,'m.'); 
hold on; h2=plot([0:511]*44.1/1024,LTqk,'r-.'); 
xlabel('frequency(kHz)');ylabel('SPL(dB)') 
legend('tonal components','thrshold in quiet') 
 
		%%%%%%%%%%%(3)NonTonalPressure%%%%%%%% 
  	for critband=1:27	%no. 0 bands don't consider 
      averagek=1;k0=0;	tmp=0; 
		  for k=boundk(critband+1):boundk(critband+2)-1 
       			averagek=averagek.*k; 
       			tmp=tmp+Xabspow(k+1); 
       			k0=k0+1; 
      end 
      if tmp~=0 
         averagek=round(averagek.^(1/k0)); 
         difference=abs( k-averagek ); 
         k0=boundk(critband+1); 
			   for k=boundk(critband+1)+1:boundk(critband+2)-1 
      			if abs(k-averagek)<difference 
        				difference=abs(k-averagek); 
        				k0=k; 
      			end 
         end 
         if tonal_flag(k0+1)==0 
			      Xnm(k0+1)=10*log10(tmp); 
			      tonal_flag(k0+1)=1; 
			   end 
      end 
   end 
figure; stem([0:511]*44.1/1024,Xnm,'.');	%stem([0:511],tonal_flag*10,'r.'); 
hold on; plot([0:511]*44.1/1024,LTqk,'r-.'); 
xlabel('frequency(kHz)');ylabel('SPL(dB)') 
legend('nontonal components','thrshold in quiet') 
 
		%%%%%%%%(4)compare with LTq() 
	for k=0:511	 
	   	if  tonal_flag(k+1)==1 & (Xnm(k+1)<LTqk(k+1))  
        		tonal_flag(k+1)=0; 
        		Xnm(k+1)=0; 
     	elseif tonal_flag(k+1)==2 & (Xtm(k+1)<LTqk(k+1))  
        		tonal_flag(k+1)=0; 
        		Xtm(k+1)=0; 
   		end 
	end 
		%%%%%%%%%(5)CompressInhalfBark() 
	for k=0:464	%band 0 don't consider 
		if tonal_flag(k+1)==2  
			j=1;	Xtmmax=Xtm;	k0=k; 
			while (k+j<464) & (CBRzk(k+j+1)-CBRzk(k+1))<=0.5 
				if tonal_flag(k+j+1)==2  
					if Xtm(k+j+1)>Xtmmax 
						tonal_flag(k0+1)=0; 
   					Xtm(k0+1)=0;  
   					Xtmmax=Xtm(k+j+1);	 
   					k0=k+j; 
   				else		 
						tonal_flag(k+j+1)=0; 
   					Xtm(k+j+1)=0; 
   				end	 
    	 end 
   		 j=j+1; 
     end 
    end 
  end 
figure; stem([0:511]*44.1/1024,Xtm);	hold on; stem([0:511]*44.1/1024,Xnm,'.'); 
xlabel('frequency(kHz)');ylabel('SPL(dB)') 
legend('tonal components','nontonal components') 
 	 
 		%%%%%%%(6)IndividMaskThreshold()	用抽选后的序号INDEX,同一i处受到的掩蔽效应加起来 
 	%%%%asign k to closest index 
 	tonal_flagindex=zeros(1,131);	Xtmindex=zeros(1,131);	Xnmindex=zeros(1,131); 
	for k=0:463 
		if tonal_flag(k+1)~=0 & tonal_flagindex(k_index(k+1)+1)~=2 
				tonal_flagindex(k_index(k+1)+1)=tonal_flag(k+1); 
				Xtmindex(k_index(k+1)+1)=Xtm(k+1); 
				Xnmindex(k_index(k+1)+1)=Xnm(k+1); 
			else  
		end 
%stem(index_k(k_index(k+1)+1),tonal_flagindex(k_index(k+1)+1)*10,'.'); 
%stem(index_k(k_index(k+1)+1),Xtmindex(k_index(k+1)+1),'c.'); 
%stem(index_k(k_index(k+1)+1),Xnmindex(k_index(k+1)+1),'g.'); 
	end 
	%%%%%%% 
  	for k=0:130 
        	LTtmsum(k+1)=0; 
        	LTnmsum(k+1)=0; 
        end 
  	for j=1:130 
        	if  tonal_flagindex(j+1)~=0 
        		if tonal_flagindex(j+1)==2 
				avtm=-1.525-0.275*CBRz(j+1)-4.5;          
			else 
				avnm=-1.525-0.175*CBRz(j+1)-0.5; 
			end 
			for k=1:130 
			  	dz=CBRz(k+1)-CBRz(j+1); 
			  	if dz>=-3 & dz<8 
		  			if -3<=dz & dz<-1  
						vf=17*(dz+1)-(0.4*X(index_k(j+1)+1)+6);    % {[-78.4,32.4)} 
	  				elseif  -1<=dz & dz<0  
						vf=(0.4*X(index_k(j+1)+1)+6)*dz;           % {[-44.4,32.4]} 
	  				elseif  0<=dz & dz<1  
						vf=-17*dz;                             % {(-17,0]} 
	  				elseif  1<=dz & dz<8  
						vf=-(dz-1)*(17-0.15*X(index_k(j+1)+1))-17; % {[-236.8,536)} 
					end 
					if tonal_flagindex(j+1)==2 
		  				LTtm=Xtmindex(j+1)+avtm+vf; 	%DB 
	  					LTtmsum(k+1)=LTtmsum(k+1)+10.^(0.1*LTtm); %ENERGY 
		  			else 
		  				LTnm=Xnmindex(j+1)+avnm+vf; 	%%%%%%%%%%%%%############## 
				  		LTnmsum(k+1)=LTnmsum(k+1)+ 10.^(0.1*LTnm); 
					end 
				end 
			end 
		end 
	end 
	 
		%%%%%%%%%%(7) global masking threshold,只计算512点  
	LTgk=zeros(1,512); 
	LTgindex=zeros(1,131); 
  	for k=1:130 
	        LTgindex(k+1)=10*log10( 10.^(0.1*LTqindex(k+1))+LTtmsum(k+1)+LTnmsum(k+1));%-normal; 
	end 
	for k=1:464 
		if k_index(k+1)~=0 
	        	LTgk(k+1)=LTgindex(k_index(k+1)+1); 
	        else 
	        	LTgk(k+1)=LTgk(k); 
	        end 
	end 
	for k=465:511 
		LTgk(k+1)=LTgk(464+1); 
	end 
	LTgk(1)=LTgk(464+1); 
figure;plot([0:511]*44.1/1024,LTgk,'r--');	hold on 
	for sb=0:31 
		LTgmin=min(LTgk(16*sb+1:16*(sb+1))); 
		for k=16*sb:16*(sb+1)-1 
			LTgk(k+1)=LTgmin; 
		end 
	end 
hold on; plot([0:511]*44.1/1024,LTgk,'k');%figure6 
xlabel('frequency(kHz)');ylabel('SPL(dB)') 
legend('global masking threshold','minimum masking threshold') 
 
 
%反归一化的PSD和掩蔽阈值 
figure;	  figure;	plot([0:511]*44.1/1024,X(1:512)-normal,'r--');	xlabel('频率(kHz)'); ylabel('SPL(dB)');hold on;				 
hold on; plot([0:511]*44.1/1024,LTgk-normal,'k');%figure8 
xlabel('frequency(kHz)');ylabel('SPL(dB)') 
legend('PSD','masking threshold') 
 
%%%%%%%%%%心理声学模型结束,输出掩蔽阈值曲线LTgk(1:512),单位分贝,后面512点与之对称,PN的FFT应在每个频率处低于该曲线 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
%%%%%%%%%%%%%%%%%%对时域整形后的PN 做FFT, 频域整形, 再IFFT%%%%%%%% 
 
 
	W1=fft(w,1024);			  
	W=20*log10(abs(W1)); 
	W(1:512)=W(1:512)+(LTgk(1:512)-normal+60.206)'; 
	%W(1)和W(513)独立不与任何对称,且W(513)没有计算 
	W(513)=W(512); 
	for k=513:1023   
		W(k+1)=W(1024-k+1); 
	end 
%plot([0:1023]/1024*44.100,W,'b..'); 
	%再向下平移直到低于掩蔽曲线 
	tmp=max(abs(W(2:512))-(LTgk(2:512)-normal+60.206)'); 
	W=W-tmp; 
  %IFFT得到最终的水印的时域 
	AI=sqrt(-1);				%%%%注意i 
	W=(10.^(W/20)).*(exp(AI*angle(W1)));  %相位采用原相位 
	w=ifft(W,1024); 
 
 
 
%学位论文画图用程序 
%(1)查表得到的44100Hz的静态掩蔽阈值图 
stem(index_k(1:131)*44.100/1024,LTqindex,'.') 
xlabel('frequency(kHz)');ylabel('SPL(dB)') 
%(2)查表得到的临界带范围 
figure;	hold on 
for k=0:28 
plot([boundk(k+1) boundk(k+1)]/1024*44.1,[0 1]) 
end 
xlabel('frequency(kHz)') 
%(3)子带的划分 
sbboundk=[0,15.5,31.5,47.5,63.5,79.5,95.5,111.5,127.5,143.5,159.5,175.5,191.5,207.5,223.5,239.5,255.5,271.5,287.5,303.5,319.5,335.5,351.5,367.5,383.5,399.5,415.5,431.5,447.5,463.5,479.5,495.5,511.5]; 
figure;	hold on 
for k=0:32 
plot([sbboundk(k+1) sbboundk(k+1)]/1024*44.1,[0 1]) 
end 
xlabel('frequency(kHz)')