psd_ex.m

用dsp解压mp3程序的算法

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% psd_ex.m - This program computes 
%            (a) periodogram
%            (b) averaged periodogram (Bartlett)
%            (c) averaged periodogram (Welch)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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)));
FFTX = fft(x1,NFFT);

% Calculate number of unique points

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

% Freq axis vector

f = (0:No_unique_pt-1)*Fs/NFFT;

% Magnitude square

Pxx = (MX.^2);

figure(2);

% Plot magnitude and phase

subplot(3,1,1);
plot(f,10*log10(Pxx));xlabel('Freq (Hz)'); ylabel('Power in dB');
axis([0 5000 -100 50])
title('The Periodogram'); grid on;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% (b) FFT for averaged peridogram - Bartlett
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

NFFT1 = 256;
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;
end    

% Averaging

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

figure(2);

% Plot magnitude and phase spectra

subplot(312);
plot(f1,10*log10(Pxx_av_bartlett));xlabel('Freq (Hz)'); ylabel('Power in dB');
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;

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;
end    

% Compute power in window function

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

% Averaging

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

figure(2);

% Plot magnitude and phase

subplot(313);
plot(f1,10*log10(Pxx_av_welch));xlabel('Freq (Hz)'); ylabel('Power in dB');
axis([0 5000 -100 50])
title('The Averaged Periodogram-Welch');grid on;