example8_3.m

用dsp解压mp3程序的算法

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%example8_3.m - FFT of sinewave in double precision,
%               single precision, and Q.15 formats 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

close all, clear all;

% Compute the FFT using double precision

M= 128;
x =sin(2*pi*(0:127)*0.1);
X = fft(x);
magX=abs(X);
phaseX=phase(X);
figure(1);
plot(magX); title('Magnitude plot using double precision');
figure(2);
plot(phaseX); title('Phase plot using double precision');


% Compute FFT using Q.15 format

Fq15=qfft('length',M,'radix',2);
Fq15.CoefficientFormat = quantizer('fixed', 'round', 'saturate',[16 15]);
Fq15.InputFormat = quantizer('fixed', 'round', 'saturate',[16 15]);
Fq15.MultiplicandFormat = quantizer('fixed', 'round', 'saturate',[16 15]);
Fq15.OutputFormat = quantizer('fixed', 'round', 'saturate',[16 15]);
Fq15.ProductFormat = quantizer('fixed', 'round', 'saturate',[32 31]);
Fq15.SumFormat = quantizer('fixed', 'round', 'saturate',[32 31]);
Fq15.ScaleValues = 1;

% Fq15.ScaleValues = [1 0.5 0.25 0.125 0.25 0.5 1];

Xq15 = fft(Fq15,x);
magXq15=abs(Xq15);
phaseXq15=phase(Xq15);
figure(1), 
plot(1:M,magX,'b',1:M,magXq15,'r'); title('Magnitude plot using double precision vs. Q.15');
legend('double precision', 'Q.15');
figure(2), 
plot(1:M, phaseX,'b',1:M,phaseXq15,'r'); title('Phase plot using double precision vs. Q.15');
legend('double precision', 'Q.15');
sum_sq_error_mag_q15=sum((magX-magXq15).^2)
sum_sq_error_phase_q15=sum((phaseX-phaseXq15).^2)

% report(Fq15)
%                        Max            Min     NOverflows    NUnderflows    NOperations
% Coefficient              1             -1              7              6            254
%       Input         0.9511        -0.9511              0              7            128
%      Output              1             -1              0              0            256
% Multiplicand             1             -1              0              9           3584
%     Product              1             -1             12              0           1792
%         Sum              2             -2            399              0           4480

% Compute FFT using Q.31 format

Fq31=qfft('length',M,'radix',2);
Fq31.CoefficientFormat = [32 31];
Fq31.InputFormat = [32 31];
Fq31.MultiplicandFormat = [32 31];
Fq31.OutputFormat = [32 31];
Fq31.ProductFormat = [64 63];
Fq31.SumFormat = [64 63];
Fq31.ScaleValues = 1;

% Fq15.ScaleValues = [1 0.5 0.25 0.125 0.25 0.5 1];

Xq31 = fft(Fq31,x);
magXq31=abs(Xq31);
phaseXq31=phase(Xq31);
figure(3),
plot(1:M,magX,'b',1:M,magXq31,'g'); title('Magnitude plot using double precision vs. Q.31');
legend('double precision', 'Q.31');
figure(4),
plot(1:M,phaseX,'b', 1:M,phaseXq31,'g'); title('Phase plot using double precision vs. Q.31');
legend('double precision', 'Q.31');

sum_sq_error_mag_q31=sum((magX-magXq31).^2)
sum_sq_error_phase_q31=sum((phaseX-phaseXq31).^2)

%qreport(Fq31)
%                        Max            Min     NOverflows    NUnderflows    NOperations
% Coefficient              1             -1              7              6            254
%       Input         0.9511        -0.9511              0              7            128
%      Output              1             -1             32              0            256
%Multiplicand              1             -1            146              0           3584
%     Product              1             -1              0              0           1792
%         Sum              2             -2            399              0           4480

    
% Compute FFT using single precion IEEE-754 format

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

Xsingle = fft(Fsingle,x);
magXsingle=abs(Xsingle);
phaseXsingle=phase(Xsingle);
figure(5), 
plot(1:M,magX,'b',1:M,magXsingle,':c'); title('Magnitude plot using double precision vs. single precision');
legend('double precision', 'single precision');
figure(6), 
plot(1:M,phaseX,'b',1:M,phaseXsingle,':c'); title('Phase plot using double precision vs. single precision');
legend('double precision', 'single precision');
qreport(Fsingle);

sum_sq_error_mag_single=sum((magX-magXsingle).^2)
sum_sq_error_phase_single=sum((phaseX-phaseXsingle).^2)

%                   Max            Min     NOverflows    NUnderflows    NOperations
% Coefficient              1             -1              0              0            254
%       Input         0.9511        -0.9511              0              0            128
%      Output          48.76         -48.76              0              0            256
% Multiplicand         30.03         -30.03              0              0           3584
%     Product          17.89         -24.12              0              0           1792
%         Sum          48.76         -48.76              0              0           4480