Matlab代码:
clear,clc,close all
L=128;dM=4;dN=8;k=0:L-1;
sigma=dM*L/(dN*2*pi);
h=(pi*(sigma))^(-0.25)*exp(-(k-0.5*(L-1)).^2/(2*sigma));
Ls=768;T=0.001;
t=(0:Ls-1)*T;
s=[sin(200*pi*t(1:256)),sin(100*pi*t(257:512)),sin(400*pi*t(513:768))];
L0=L+Ls-dM;s1=[zeros(1,L-dM),s];h1=[h,zeros(1,Ls-dM)]/norm(h);
h=[h1,h1];s=[s1,s1];N=L0/dN;M=L0/dM;H=zeros(dM*dN,L0);
for m=0:dN-1
for n=0:dM-1
for k=0:L0-1
H(n+m*dM+1,k+1)=h(floor(k+m*N+1))*exp(-2*pi*i*n*k/dM);
end
end
end
mu=[dM/N,zeros(1,dM*dN-1)]';
rg=H'*inv(H*H')*mu;
g=rg'/norm(rg);g=g(1:L);
figure(1),
plot((0:L-1),h1(1:L)),
hold on
plot((0:L-1),real(g),'r--'),
xlabel('k'),ylabel('amplitude'),title('标架与对偶标架')
legend('Gabor标架','Gabor对偶标架')
Cmn=zeros(M,L0);
for m=0:M-1
for n=0:L0-1
sg(n+1)=s(n+m*dM+1)*rg(n+1);
end
Cmn(m+1,:)=fft(sg);
end
isgabor=zeros(M,L0);
for i=1:M
isgabor(i,:)=ifft(Cmn(i,:));
for j=1:L0
isg(i,(i-1)*dM+j)=h1(j)*isgabor(i,j);
end
end
isg=ones(1,M)*isg*N;is=isg(1:L0);error=abs(s1)-abs(is);figure(2);subplot(221),
plot(t,s1(L-dM+1:L0)),
xlabel('time/s'),ylabel('amplitude'),title('input signal');
axis tight;
subplot(223),
plot(t,real(is(L-dM+1:L0))),
xlabel('time/s'),ylabel('amplitude'),title('construct signal');
axis tight;
subplot(224),
plot(t,error(L-dM+1:L0));
xlabel('time/s'),ylabel('amplitude'),title('误差');
axis tight;
spectrogram=abs(Cmn((L-dM)/dM:(Ls-dM)/dM,:)');
maxspec=max(max(spectrogram));
minspec=min(min(spectrogram));
colormap(1-gray(256)),
subplot(222),
image(t,(0:L0/2-1)/(L0*T),256*(spectrogram(1:L0/2,:)-minspec)/(maxspec-minspec));
xlabel('time/ms'),ylabel('frequency/Hz'),title('时频域')
axis('xy'),colorbar
clear
L=128;
dN=8;
dM=16;
M=L/dM;N=L/dN;
k=0:L-1;
seigma2=dM*L/(dN*2*pi);
g=(pi*seigma2).^(-1/4)*exp(-(k-(L-1)/2).^2/(2*seigma2));
h=zeros(1,L)';
mu=[dM/N zeros(1,dM*dN-1)]';
H=zeros(dM*dN,L);
for m=0:dN-1
for n=0:dM-1
for k=0:L-1
H(n+m*dM+1,k+1)=g(floor(mod(k+m*N,L))+1)*exp(-2*pi*i*n*k/dM);
end
end
end
h=H'*inv((H)*H')*mu;
h=h/norm(h);
figure(3)
plot([0:L-1],g),hold on
plot([0:L-1],real(h),'r--')
xlabel('L=256,dM=10,dN=12,紧标架'),ylabel('Amplitude'),title('标架与对偶标架');
legend('标架',' 对偶标架')