% Copyleft (C) Forrest Sheng Bao 1984-2009
% with Dept. of Electrical and Computer Engineering
% and  Dept. of Computer Science
% at Texas Tech University, Lubbock, Texas, USA
% http://fsbao.net
%
% This program shows how Welch method can change the spectrum
%
% This program is a free software.
% It is licensed in GNU GENERAL PUBLIC LICENSE (GPL) v3 or later.
% So this program comes with ABSOLUTELY NO WARRANTY.
%
% You have the rights to freely COPY, MODIFY and REDISTRIBUTE this program
% under the terms of GNU GPL v3 or later, including KEEP the modified
% and/or redistributed code in GNU GPL v3 or later.
%
% If you do not know what GNU GPL is, please visit
% www.gnu.org/copyleft/gpl.html

clear all;
close all;

Fs = 1000;%sampling frequency

y = sig_noise([140 180], [-10 -10], 256);% generating the signal
L = length(y);%length of the signal


figure();
plot(1:L,y);%plot the signal in time domain
axis tight
axis 'auto y'
title('256-point signal containing two sines@140 and 180 Hz, both with SNR of -10db ')
xlabel('time (milliseconds)')
ylabel('amplitude')

% Below is direct FFT and Welch method

NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);%only one side of the spectrum

figure();

% Plot single-sided amplitude spectrum.
plot(f,2*abs(Y(1:NFFT/2+1)))
title('Single-sided amplitude spectrum of generated signal')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')

for overlap_percent = [0 0.5 0.99]
    NFFT = 64; % FFT points
    SegLen = 64; % length of each segment
    overlap = floor(SegLen*overlap_percent); % overlap percentage

    [PS,f] = pwelch(y, hamming(SegLen), overlap, NFFT, Fs); % doing Welch method
    figure();
    plot(f,PS)
    title([num2str(SegLen) '-sample Hamming window, ' num2str(NFFT) '-point FFT, ' num2str(overlap_percent) '% overlap'])
    xlabel('Frequency (Hz)')
    ylabel('Welch method PSD, |Y(f)|')
end

Published with MATLAB® 7.8