psd_ex_quant.m

用dsp解压mp3程序的算法

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% psd_ex_quant.m - This m-file using different formats to compute 
%                     (a) periodogram
%                     (b) averaged periodogram (Bartlett)
%                     (c) averaged periodogram (Welch)
%                  in Q.15 implementation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all; close all;

% Generate and mix 3 sinewaves

F1 = 1000; F2 = 1200; F3 = 1400;
Fs = 10000;
T = 0:1:499;
x1 = sin(2*pi*T*F1/Fs) + 0.5*sin(2*pi*T*F2/Fs) + 0.25*sin(2*pi*T*F3/Fs);
figure(1);plot(x1); title('Sine waves');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% (a) FFT for averaged peridogram - Periodogram
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

NFFT = 2^(nextpow2(length(x1)));

% Compute FFT using Q.15 format

Fscale=qfft('length',NFFT,'radix',2);
Fscale.CoefficientFormat = [16 15];
Fscale.InputFormat = [16 15];
Fscale.MultiplicandFormat = [16 15];
Fscale.OutputFormat = [16 15];
Fscale.ProductFormat = [32 31];
Fscale.SumFormat = [32 31];
Fscale.ScaleValues = [0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5];

% Single precision floating-point works better

Fsingle=qfft('length',NFFT,'radix',2);
Fsingle.quantizer = {'single'};

FFTXq = fft(Fscale,x1);
FFTX = fft(x1, NFFT);

% Calculate number of unique points

No_unique_pt = ceil((NFFT+1)/2);
FFTX = FFTX(1:No_unique_pt);
FFTXq = FFTXq(1:No_unique_pt);
MX = abs(FFTX)/NFFT;
MXq = abs(FFTXq);

f = (0:No_unique_pt-1)*Fs/NFFT;  % frequency axis vector
Pxx = (MX.^2);                   % magnitude square
Pxxq = (MXq.^2);

figure(2);                       % Plot magnitude and phase
plot(f,10*log10(Pxx),'r',f,10*log10(Pxxq),'g');xlabel('Freq (Hz)'); ylabel('Power in dB');
legend('double precision','Q.15');
axis([0 5000 -100 50]); xlabel('Freq (Hz)'); ylabel('Power in dB');
title('The Periodogram');grid on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% (b) FFT for averaged peridogram - Bartlett
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT1 = 256;

% Compute FFT using Q.15 format

Fscale=qfft('length',NFFT1,'radix',2);
Fscale.CoefficientFormat = [16 15];
Fscale.InputFormat = [16 15];
Fscale.MultiplicandFormat = [16 15];
Fscale.OutputFormat = [16 15];
Fscale.ProductFormat = [32 31];
Fscale.SumFormat = [32 31];
Fscale.ScaleValues = [0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5];

% Single precision floating-point works better
Fsingle=qfft('length',NFFT1,'radix',2);
Fsingle.quantizer = {'single'};

K = round(length(x1)/NFFT1);

% Compute segment FFT

for seg = 0:1:K-1
    if NFFT1+(seg*NFFT1) > length(x1)
        x_seg = x1(1+(seg*NFFT1):1:length(x1));
    else    
        x_seg = x1(1+(seg*NFFT1):1:NFFT1+(seg*NFFT1));
    end    
    X(seg+1,:) = (abs(fft(x_seg,NFFT1))/NFFT1).^2;
    Xq(seg+1,:) = (abs(fft(Fscale,x_seg))).^2;
end    

% Averaging

X_add = zeros(1,NFFT1);
X_addq = zeros(1,NFFT1);
for seg = 1:K
    X_add = X(seg,:)+X_add;
    X_addq = Xq(seg,:)+X_addq;
end    
Pxx_av_bartlett = X_add./K;
Pxx_av_bartlettq = X_addq./K;
No_unique_pt1 = ceil((NFFT1+1)/2);
Pxx_av_bartlett = Pxx_av_bartlett(1,1:No_unique_pt1);
Pxx_av_bartlettq = Pxx_av_bartlettq(1,1:No_unique_pt1);
f1 = (0:No_unique_pt1-1)*Fs/NFFT1;

figure(3);             % Plot magnitude and phase
plot(f1,10*log10(Pxx_av_bartlett),'r',f1,10*log10(Pxx_av_bartlettq),'g');xlabel('Freq (Hz)'); ylabel('Power in dB'); hold on;
legend('double precision','Q.15');
axis([0 5000 -100 50])
title('The Averaged Periodogram-Bartlett'); grid on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% (c) FFT for averaged peridogram - Welch
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT1 = 256;
overlap = NFFT1/2;          % 50% overlapping
KW = ceil((length(x1)-NFFT1)/overlap+1);

% Compute segment FFT

x_seg1 = x1(1:1:NFFT1).*hamming(NFFT1)';
X1(1,:) = (abs(fft(x_seg1,NFFT1))/NFFT1).^2;
X1q(1,:) = (abs(fft(Fscale,x_seg1))).^2;

for seg = 1:1:KW-1
    if NFFT1+(seg*NFFT1-overlap) > length(x1)
        x_seg1 = x1(1+(seg*(NFFT1-overlap)):1:length(x1));
        x_seg1 = [x_seg1 zeros(1,(NFFT1+seg*(NFFT1-overlap))-length(x1))];
        x_seg1 = x_seg1.*hamming(NFFT1)';
    else    
        x_seg1 = x1(1+(seg*(NFFT1-overlap)):1:NFFT1+(seg*(NFFT1-overlap))).*hamming(NFFT1)';
    end    
    X1(seg+1,:) = (abs(fft(x_seg1,NFFT1))/NFFT1).^2;
    X1q(seg+1,:) = (abs(fft(Fscale,x_seg1))).^2;
end    

% Compute power in window function

U = sum((hamming(NFFT1)).^2)/NFFT1;

% Averaging
X_add1 = zeros(1,NFFT1);
X_add1q = zeros(1,NFFT1);
for seg = 1:KW
    X_add1 = X1(seg,:)+X_add1;
    X_add1q = X1q(seg,:)+X_add1q;
end    
Pxx_av_welch = X_add1./(U*KW);
Pxx_av_welchq = X_add1q./(U*KW);
No_unique_pt1 = ceil((NFFT1+1)/2);
Pxx_av_welch = Pxx_av_welch(1,1:No_unique_pt1);
Pxx_av_welchq = Pxx_av_welchq(1,1:No_unique_pt1);
f1 = (0:No_unique_pt1-1)*Fs/NFFT1;

figure(4);                % Plot magnitude and phase

plot(f1,10*log10(Pxx_av_welch),'r',f1,10*log10(Pxx_av_welchq),'g');xlabel('Freq (Hz)'); ylabel('Power in dB'); hold on;
legend('double precision','Q.15');
axis([0 5000 -100 50])
title('The Averaged Periodogram-Welch'); grid on;