108 lines
2.4 KiB
Mathematica
108 lines
2.4 KiB
Mathematica
|
%% Clear Workspace and Close figures
|
||
|
|
clear; close all; clc;
|
||
|
|
|
||
|
|
%% Intialize Laplace variable
|
||
|
|
s = zpk('s');
|
||
|
|
|
||
|
|
% Signal's PSD
|
||
|
|
% We load the PSD of the signal we wish to replicate.
|
||
|
|
|
||
|
|
load('./mat/dist_psd.mat', 'dist_f');
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
% We remove the first value with very high PSD.
|
||
|
|
|
||
|
|
dist_f.f = dist_f.f(3:end);
|
||
|
|
dist_f.psd_gm = dist_f.psd_gm(3:end);
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
% The PSD of the signal is shown on figure ref:fig:psd_original.
|
||
|
|
|
||
|
|
figure;
|
||
|
|
hold on;
|
||
|
|
plot(dist_f.f, dist_f.psd_gm)
|
||
|
|
hold off;
|
||
|
|
xlabel('Frequency [Hz]');
|
||
|
|
ylabel('Power Spectral Density');
|
||
|
|
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
|
||
|
|
xlim([0.1, 500]);
|
||
|
|
|
||
|
|
% Algorithm
|
||
|
|
% We define some parameters that will be used in the algorithm.
|
||
|
|
|
||
|
|
Fs = 2*dist_f.f(end); % Sampling Frequency of data is twice the maximum frequency of the PSD vector [Hz]
|
||
|
|
N = 2*length(dist_f.f); % Number of Samples match the one of the wanted PSD
|
||
|
|
T0 = N/Fs; % Signal Duration [s]
|
||
|
|
df = 1/T0; % Frequency resolution of the DFT [Hz]
|
||
|
|
% Also equal to (dist_f.f(2)-dist_f.f(1))
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
% We then specify the wanted PSD.
|
||
|
|
|
||
|
|
phi = dist_f.psd_gm;
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
% We create amplitudes corresponding to wanted PSD.
|
||
|
|
|
||
|
|
C = zeros(N/2,1);
|
||
|
|
for i = 1:N/2
|
||
|
|
C(i) = sqrt(phi(i)*df);
|
||
|
|
end
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
% Finally, we add some random phase to =C=.
|
||
|
|
|
||
|
|
theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
|
||
|
|
|
||
|
|
Cx = [0 ; C.*complex(cos(theta),sin(theta))];
|
||
|
|
Cx = [Cx; flipud(conj(Cx(2:end)))];;
|
||
|
|
|
||
|
|
% Obtained Time Domain Signal
|
||
|
|
% The time domain data is generated by an inverse FFT.
|
||
|
|
|
||
|
|
% The =ifft= Matlab does not take into account the sampling frequency, thus we need to normalize the signal.
|
||
|
|
|
||
|
|
u = N/sqrt(2)*ifft(Cx); % Normalisation of the IFFT
|
||
|
|
t = linspace(0, T0, N+1); % Time Vector [s]
|
||
|
|
|
||
|
|
figure;
|
||
|
|
plot(t, u)
|
||
|
|
xlabel('Time [s]');
|
||
|
|
ylabel('Amplitude');
|
||
|
|
xlim([t(1), t(end)]);
|
||
|
|
|
||
|
|
% PSD Comparison
|
||
|
|
% We duplicate the time domain signal to have a longer signal and thus a more precise PSD result.
|
||
|
|
|
||
|
|
u_rep = repmat(u, 10, 1);
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
% We compute the PSD of the obtained signal with the following commands.
|
||
|
|
|
||
|
|
nx = length(u_rep);
|
||
|
|
na = 16;
|
||
|
|
win = hanning(floor(nx/na));
|
||
|
|
|
||
|
|
[pxx, f] = pwelch(u_rep, win, 0, [], Fs);
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
% Finally, we compare the PSD of the original signal and the obtained signal on figure ref:fig:psd_comparison.
|
||
|
|
|
||
|
|
figure;
|
||
|
|
hold on;
|
||
|
|
plot(dist_f.f, dist_f.psd_gm, 'DisplayName', 'Original PSD')
|
||
|
|
plot(f, pxx, 'DisplayName', 'Computed')
|
||
|
|
hold off;
|
||
|
|
xlabel('Frequency [Hz]');
|
||
|
|
ylabel('Power Spectral Density');
|
||
|
|
set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
|
||
|
|
legend('location', 'northeast');
|
||
|
|
xlim([0.1, 500]);
|