Math Forum Discussions

Robert Casey

Posts: 4
Registered: 12/7/04

short m-file for waveform -> freq display jfrz.m 1/1
Posted: Aug 5, 1996 4:29 PM

Plain Text

A small m-file we wrote up to display the spectra of signals of digital TV work,
also works for other kinds of digitized analog signals.

function jfrz(n,d,f,mode,gain,color);
% this M file computes freq response for the input
% n = the input numerator
% d = the input denominator
% f = sampling frequency used calculate x scale reference
% mode = picture view mode
% 0 -- erase screen and draw new waveform
% 1 -- add to existing plot
% gain = how dc gain is implimented
% 0 -- normalize max value to 0 db
% 1 -- no normalizing
% color = the view color
% 0 -- black background yellow data
% 1 -- white background black data
%function jfrz(n,d,f,mode,gain,color);
%if (nargin == 1) d=1; f=27e6; mode=0; gain=0; color=0; end
%if (nargin == 2) f=27e6; mode=0; gain=0; color=0; end
%if (nargin == 3) mode=0; gain=0; color=0; end
%if (nargin == 4) gain=0; color=0; end
%if (nargin == 5) color=0; end

npoints = 4096;
min=10e-8;

if (nargin == 1) d=1; f=27e6; mode=0; gain=0; color=0; end
if (nargin == 2) f=27e6; mode=0; gain=0; color=0; end
if (nargin == 3) mode=0; gain=0; color=0; end
if (nargin == 4) gain=0; color=0; end
if (nargin == 5) color=0; end

[h,w]= freqz(n,d,npoints);
a = abs(h);
if (gain == 0)
a=a/max(a);
end;
tt = find(a < min);
a(tt)=min*ones(1,size(tt,1));
a=20*log10(a);

%fsc = round(2*npoints*3.57954e6/f);
%a(fsc)=0;

f=f/2;
fr = [0:f/(npoints-1):f]';

if (mode == 0)
clf reset
axes('position',[0.06 0.10 0.9 0.85]);
end;
if (mode == 1)
hold on;
end;
if (color == 0)
set(gcf,'color',[0 0 0]);
plot(fr,a,'y');
set(gca, 'color',[0 0 0]);
set(gca,'Xcolor',[1 1 1]);
set(gca,'Ycolor',[1 1 1]);
end;
if (color == 1)
set(gcf,'color',[1 1 1]);
plot(fr,a,'k');
set(gca, 'color',[1 1 1]);
set(gca,'Xcolor',[0 0 1]);
set(gca,'Ycolor',[0 0 1]);
end;
grid

you'll also need freqz.m (from signal toolbox, below for reference only):

function [h,w] = freqz(b,a,n,dum)
%FREQZ Z-transform digital filter frequency response. When N is an integer,
% [H,W] = FREQZ(B,A,N) returns the N-point frequency vector W and the
% N-point complex frequency response vector G of the filter B/A:
% -1 -nb
% jw B(z) b(1) + b(2)z + .... + b(nb+1)z
% H(e) = ---- = ----------------------------
% -1 -na
% A(z) 1 + a(2)z + .... + a(na+1)z
% given numerator and denominator coefficients in vectors B and A. The
% frequency response is evaluated at N points equally spaced around the
% upper half of the unit circle. To plot magnitude and phase of a filter:
% [h,w] = freqz(b,a,n);
% mag = abs(h); phase = angle(h);
% semilogy(w,mag), plot(w,phase)
% FREQZ(B,A,N,'whole') uses N points around the whole unit circle.
% FREQZ(B,A,W) returns the frequency response at frequencies designated
% in vector W, normally between 0 and pi. (See LOGSPACE to generate W).
% See also YULEWALK, FILTER, FFT, INVFREQZ, and FREQS.

% J.N. Little 6-26-86
% Revised 6-7-88 JNL, 9-12-88 LS
% Copyright (c) 1985-1989 by the MathWorks, Inc.

a = a(:).';
b = b(:).';
na = max(size(a));
nb = max(size(b));
nn = max(size(n));
s = 2;
if nargin == 4
s = 1;
end
if nn == 1
w = (2*pi/s*(0:n-1)/n)';
else
w = n;
n = nn;
end
if (nn == 1)
h = (fft([b zeros(1,s*n-nb)]) ./ fft([a zeros(1,s*n-na)])).';
h = h(1:n);
else
% Frequency vector specified. Use Horner's method of polynomial
% evaluation at the frequency points and divide the numerator
% by the denominator.
a = [a zeros(1,nb-na)]; % Make sure a and b have the same length
b = [b zeros(1,na-nb)];
s = exp(sqrt(-1)*w);
h = polyval(b,s) ./ polyval(a,s);
end