c8_psdexample.m

经典书籍《通信系统仿真原理与无线应用》chap8部分课后题仿真源代码

% File: c8_PSDexample.m
settle = 100;                  % ignore transient
fs = 1000;                    % sampling frequency
N = 50000;                   % size of data record
f = (0:(N-1))*fs/N;              % frequency scale
[b,a] = cheby1(5,5,0.1);        % filter
NN = N+settle;                % allow transient to die
in = randn(1,NN);             % random input
out = filter(b,a,in);             % filter output
out = out((settle+1):NN);       % strip off initial samples
window = hanning(N)';         % set window function
winout = out.*window;         % windowed filter output
fout = abs(fft(winout,N)).^2;     % transform and square mag
U = sum(window.*window);      % window energy
flout = fout/U;               % scale spectrum
psdl = 10*log10(abs(flout));  % log scale
subplot(2,1,1)
plot(f(1:5000),psdl(1:5000))
grid; axis([0 100 -70 10]);
xlabel('Frequency, Hz')
ylabel('PSD')
%
K = 25;                   % number of segments
M = N/K;                  % block size
fK = (0:(M-1))*fs/M;      % frequency scale
d = zeros(1,M);           % initialize vector
psdk = zeros(1,M);        % initialize vector
%window = hanning(M)';     % set window function
window = boxcar(M)';
U = sum(window.*window);  % window energy
for k=1:K
for j=1:M
index = (k-1)*M+j;
d(j) = out(index);
end
dwin = d.*window;
psdk = (abs(fft(dwin,M)).^2)/U + psdk;
end
psd2 = 10*log10(psdk/K);
subplot(2,1,2)
plot(fK(1:250),psd2(1:250))
grid; axis([0 100 -70 10]);
xlabel('Frequency, Hz')
ylabel('PSD')
% End of script file.