fg_10_11.m

英文书《Digital Signal Processing with Examples in MATLAB》附带的MATLAB实例

% fg_10_11          
K=20000; T=1/40000; skp=5000;
% Copy the speech signal.
fid=fopen(['signal_3.txt']);
s=fscanf(fid,'%7d%8d%8d%8d%8d%8d%8d%8d%8d%8d\n',K);
fclose(fid);
s=round(127*s(skp+1:K)'/max(abs(s(skp+1:K))));
N=length(s);
% Plot s and spectra.
sp_fig(1);
% Signal
axes('position',[.15,.7,.7,.27]);
set(gca,'fontsize',16);
line([0:N-1]*T,s,'color','r'); grid on;
axis([0 (N-1)*T -150 150]);
xlabel('time (s)'); ylabel('s(t)');
% DFT amplitude.
axes('position',[.06,.12,.4,.43]);
set(gca,'fontsize',16);
S=fft(s);
amp=abs(S(1:N/2+1));
line([0:N/2]*.001/(N*T),amp,'color','b'); grid on;
axis([0 (N/2)*.001/(N*T) 0 40000]);
xlabel('Frequency (KHz)'); ylabel('DFT amplitude');
% DCT magnitude(twice as many points).
axes('position',[.56,.12,.4,.43]);
set(gca,'fontsize',16);
S=sp_dct(s);
line([0:N-1]*.001/(2*N*T),abs(S),'color','b'); grid on;
axis([0 (N/2)*.001/(N*T) 0 500]);
xlabel('Frequency (KHz)'); ylabel('DCT amplitude');
% Check.
Y_sds=[S(1:6); S(N-1)]'
Y=dct(s);
Y_Matlab=[Y(1:6) Y(N-1)]