This document contains code to generate a Gaussian process, filter it with a low-pass filter, downsample the filtered signal, and analyze the power and statistics of the original and downsampled signals. It generates the Gaussian process, applies a filter, convolves the output with an impulse response. It then calculates autocorrelations, power, power spectral density of the original signal. It downsamples the signal, recalculates these metrics, and estimates how power and other statistics vary with sampling period by introducing more downsampling. It finally calculates the mean square error between the original and downsampled power values.
This document contains code to generate a Gaussian process, filter it with a low-pass filter, downsample the filtered signal, and analyze the power and statistics of the original and downsampled signals. It generates the Gaussian process, applies a filter, convolves the output with an impulse response. It then calculates autocorrelations, power, power spectral density of the original signal. It downsamples the signal, recalculates these metrics, and estimates how power and other statistics vary with sampling period by introducing more downsampling. It finally calculates the mean square error between the original and downsampled power values.
This document contains code to generate a Gaussian process, filter it with a low-pass filter, downsample the filtered signal, and analyze the power and statistics of the original and downsampled signals. It generates the Gaussian process, applies a filter, convolves the output with an impulse response. It then calculates autocorrelations, power, power spectral density of the original signal. It downsamples the signal, recalculates these metrics, and estimates how power and other statistics vary with sampling period by introducing more downsampling. It finally calculates the mean square error between the original and downsampled power values.
This document contains code to generate a Gaussian process, filter it with a low-pass filter, downsample the filtered signal, and analyze the power and statistics of the original and downsampled signals. It generates the Gaussian process, applies a filter, convolves the output with an impulse response. It then calculates autocorrelations, power, power spectral density of the original signal. It downsamples the signal, recalculates these metrics, and estimates how power and other statistics vary with sampling period by introducing more downsampling. It finally calculates the mean square error between the original and downsampled power values.
Download as DOC, PDF, TXT or read online from Scribd
Download as doc, pdf, or txt
You are on page 1/ 15
Order #FO7C70BFCC81
Buyer: nicholassalis December 14, 2018
Version 1.0 (December 15, 2018)
PLEASE DELETE THIS FIRST PAGE IF THIS IS THE FINAL VERSION
Input data Generate Gaussian process Filter with 1KHz cutof Generate impulse response Generate output of the filter by convolution between the process realizations and the impulse response 1.1 Autocorrelation plot 1.2 Calculate Power 1.3 Power spectral density - TO CHOOSE METHOD 2. Sample with n*T -> basically this means decimation. Let's try by a factor of 2 3 Estimate power based on a limited number of N samples, for diferent sampling periods 4. Mean Square Error
% Clean Start clear close all
Input data Bw = 1e3; % Signal bandwidth Fs = 10*Bw;% Sampling frequency must be at least 2*Bw. Choose oversampling dt = 1/Fs; % Sampling period - useful for constructing the time array N = 2e4; % Number of samples R = 1e2; % Number of realizations m = 1; % mean value sigma = 1; % standard deviation
Generate Gaussian process
x = sigma*randn(N,R) + m; t = 0:dt:dt*(N-1); t = t-t(end)/2; t = t'; figure(1) plot(t,x) xlabel('t[s]') ylabel('x(t)') title(['All the ',num2str(R),' realizations vs time']) Filter with 1KHz cutof Wn = Bw/(Fs/2); ord = 200; [B,A] = fir1(ord,Wn,'low'); x_filtered = filtfilt(B,A,x); figure(2) plot(t,x_filtered) xlabel('t[s]') ylabel('x_{filtered}(t)') title(['All the ',num2str(R),' filtered realizations vs time']) Generate impulse response h = exp(-t/0.1).*heaviside(t); figure(3) plot(t,h) xlabel('t[s]') ylabel('h(t)') title('Filter impulse response') Generate output of the filter by convolution between the process realizations and the impulse response y = zeros(N,R); for i = 1:R y(:,i) = conv(x(:,i),h,'same'); end figure(4) plot(t,y) xlabel('t[s]') ylabel('y(t)') title(['All the ',num2str(R),' filter outputs vs time']) 1.1 Autocorrelation plot figure(5) hold on for i = 1:R [acf,lags,bounds] = autocorr(y(:,i),N-1); plot(lags,acf) end xlabel('lags(y(t))') ylabel('autocorr(y)') title(['All the ',num2str(R),' autocorrelations']) hold off 1.2 Calculate Power power_original = (1/N)*(sum(y.^2)); figure(6) plot(power_original) xlabel('realizations') ylabel('Power') title('Power of y(t) vs realizations') 1.3 Power spectral density - TO CHOOSE METHOD dF = Fs/N; F = -Fs/2:dF:Fs/2-dF; [PSD,F] = pwelch(y,[],[],F,Fs); figure(7) plot(F,PSD) xlabel('Frequency [Hz])') ylabel('PSD') title(['All the ',num2str(R),' power spectral densities']) % or Y = fft(y); Y = Y./(N); Y = Y.*conj(Y); Sy = mean(Y,2); SSy = fftshift(Y); figure(8) plot(F,SSy) xlabel('Frequency [Hz]') ylabel('PSD') title(['All the ',num2str(R),' power spectral densities']) 2. Sample with n*T -> basically this means decimation. Let's try by a factor of 2 f = 2; Fs2 = Fs/f; T = 1/Fs2; N2 = round(N/f); y2 = zeros(N2,R); figure(9) hold on for i = 1:R y2(:,i) = decimate(y(:,i),f); [acf,lags,bounds] = autocorr(y2(:,i),N2-1); plot(lags,acf) end xlabel('lags(y2(t))') ylabel('autocorr(y2)') title(['All the ',num2str(R),' autocorrelations']) hold off % Recalculate power for the resampled(decimated) signal power_decimated = (1/length(y2) )* sum(y2.^2); figure(10) plot(power_decimated) xlabel('realizations') ylabel('Resampled power') title('Resampled power of y(t) vs realizations') % Recalculate power spectral density for the resampled(decimated signal) dF = Fs2/N2; F = -Fs2/2:dF:Fs2/2-dF; [PSD,F] = pwelch(y2,[],[],F,Fs2); figure(11) plot(F,PSD) xlabel('Frequency [Hz]') ylabel('PSD') title(['All the resampled ',num2str(R),' power spectral densities']) 3 Estimate power based on a limited number of N samples, for diferent sampling periods f = [1:10]; Fs2 = Fs./f; % Sampling frequency decreases T2 = 1./Fs2; % Sampling period increases N2 = round(N./f); % total number of samples decreases Nmax = N;
for j = 1:length(f) for i = 1:R y2 = decimate(y(1:Nmax,i),f(j)); % Power powerN(j,i) = (1/length(y2))*(sum(y2.^2)); % Average value mean_N(j,i) = mean(y2); % Variance var_N(j,i) = var(y2); end end figure(12) plot(powerN') xlabel('realizations') ylabel('Power') title('Power of y(t) vs realizations') figure(13) plot(mean_N') xlabel('realizations') ylabel('Mean value') title('Mean value of y(t) vs realizations') figure(14) plot(var_N') xlabel('realizations') ylabel('Variance') title('Variance of y(t) vs realizations') 4. Mean Square Error T = f*dt; errors = (power_original - powerN).^2; mean_sq_err = sum(errors')/R; figure(15) plot(T,mean_sq_err) xlabel('Sampling period') ylabel('Mean Square Error') title('Mean Square Error vs sampling period') axis('tight')
