Skip to content

Instantly share code, notes, and snippets.

@rstemmer
Created February 23, 2013 20:27
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save rstemmer/5021206 to your computer and use it in GitHub Desktop.
Save rstemmer/5021206 to your computer and use it in GitHub Desktop.
# source "H.octave"
function [magnitude, phase] = H(omega, omegac, R)
# numerator
N = 1 - e.^(-2 * omega * j);
# denominator
D = 1 - 2 * R * cos(omegac) * e.^(-i * omega) + R^2 * e.^(-2 * omega * j);
F = N ./ D;
magnitude = abs(F);
phase = angle(F);
endfunction
fc = 2.3041E6; # center frequency
fs = 115.2E6; # sampling frequency
omegac = 2 * pi * fc / fs;
im = [-1:.03:1];
re = [-1:.03:1];
for index_im = 1:size(im,2)
for index_re = 1:size(re,2)
c = re(index_re) + i * im(index_im);
if abs(c) <= 1
[m,p] = H(angle(c), omegac, abs(c));
if m > 10
m = 10;
endif
else
m = -1;
endif
mag(index_im, index_re) = m;
endfor
endfor
hold on
axis([-1, 1, -1, 1, 0, max(max(mag))]);
mesh(re, im, mag)
t = linspace(-pi, pi, 750)';
circle_im = sin(t);
circle_re = cos(t);
c = circle_re + i * circle_im;
[m,p] = H(angle(c), omegac, 0.99);
fig = plot3(circle_re, circle_im, m, "r");
set(fig, "linewidth", 2);
fig = plot(circle_re, circle_im, "k");
set(fig, "linewidth", 2);
hold off
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment