1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
|
% impulse_noise
% David Rowe May 2017
%
% Experiments with impulsive noise and HF radio
format;
more off;
rand('seed',1)
% DFT function ------------------------------------------------
% note k is on 0..K-1 format, unlike Octave fft() which is 1..K
function H = calc_H(k, K, a, d)
L = length(d);
H = 0;
for i=1:L
H += a(i)*exp(-j*2*pi*k*d(i)/K);
end
endfunction
% -----------------------------------------
% PWM noise simulation
% -----------------------------------------
function pwm_noise
Fs = 10E6; % sample rate of simulation
Fsig = 1E6; % frequency of our wanted signal
Fpwm = 255E3; % switcher PWM frequency
T = 1; % length of simulations in seconds
Nsam = T*Fs;
Nsamplot = 200;
Apwm = 0.1;
Asig = -40; % attenuation of wanted signal in dB
% generate an impulse train with jitter to simulate switcher noise
pwm = zeros(1,Fs);
Tpwm = floor(Fs/Fpwm);
pulse_positions_pwm = Tpwm*(1:T*Fpwm) + round(rand(1,T*Fpwm));
h_pwm = zeros(1,Nsam);
h_pwm(pulse_positions_pwm) = Apwm;
h_pwm = h_pwm(1:Nsam);
% add in wanted signal and computer amplitude spectrum
s = 10^(Asig/20)*cos(2*pi*Fsig*(1:Nsam)/Fs);
h = h_pwm+s;
H = fft(h);
Hdb = 20*log10(abs(H)) - 20*log10(Nsam/2);
figure(1); clf;
subplot(211)
plot(h(1:Nsamplot));
subplot(212)
plot(Hdb(1:Nsam/2));
axis([0 T*2E6 -120 0]); xlabel('Frequency Hz'); ylabel('Amplityude dBV'); grid;
printf("pwm rms: %f signal rms: %f noise rms\n", std(h_pwm), std(s));
endfunction
% -----------------------------------------
% Single pulse noise simulation
% -----------------------------------------
function pulse_noise
% set up short pulse in wide window, consisting of two samples next
% to each other
K = 1024;
a(1) = a(2) = 1; d(1) = 10; d(2) = d(1)+1;
h = zeros(1,K);
h(d(1)) = a(1);
h(d(2)) = a(2);
% mag and phase spectrum, mag spectrum changes slowly
figure(2); clf;
Hfft = fft(h);
subplot(311)
stem(h(1:100));
axis([1 100 -0.2 1.2]);
subplot(312)
plot(abs(Hfft(1:K/2)),'+');
title('Magnitude');
subplot(313)
plot(angle(Hfft(1:K/2)),'+');
title('Phase');
% simple test to estimate H(k+1) from H(k) --------------------
% brute force calculation
k = 300;
H = zeros(1,K);
H(k-1) = calc_H(k-1, K, a, d);
H(k) = calc_H(k, K, a, d);
H(k+1) = calc_H(k+1, K, a, d);
% calculation of k+1 from k using approximation that {d(i)} are
% close together compared to M, i.e it's a narrow pulse (assumes we
% can estimate d using other means)
Hk1_ = exp(-j*2*pi*d(1)/K)*H(k);
% plot zoomed in version around k to compare
figure(3); clf;
plot(H(k-1:k+1),'b+','markersize', 10, 'linewidth', 2);
hold on; plot(Hk1_,'g+','markersize', 10, 'linewidth', 2); hold off;
title('H(k-1) .... H(k+1)');
printf("H(k+1) match: %f dB\n", 20*log10(abs(H(k+1) - Hk1_)));
endfunction
% Run various simulations here ---------------------------------------------
%pwm_noise
pulse_noise
|
