wbfm.m

很多matlab的源代码

function [y,yi] = wbfm(ty,beta,m)
% WBFM Percent Power in a Wideband FM Signal.
%
%	[Y,YC] = WBFM(TY,BETA,M) computes percent power in a wideband FM signal 
%	TY = modulation type: 'si' for sine wave, 'tr' for triangular wave, 
%	'st' for sawtooth, and 'sq' for square wave. 
%	BETA = modulation index.
%	For integer M, Y = fraction of total power in M harmonics 
%	For 0< M < 1, Y = number of harmonics containing 100M% of total power.
%	[DEFAULT: M = fix(BETA)]
%	YC is the cumulative power through M harmonics.
%
%	The spectrum is plotted if no output arguments are given.
%
% 	USAGE: The number of harmonics N that contain >75% of total power 
%	in a sawtooth FM  signal with beta= 8 is: >>N = wbfm('st',8,0.75)
%
%	WBFM (with no input arguments) invokes the following example:
%
%	% Plot the spectrum for sinusoidal FM with beta=10
%	  >>wbfm('si',10);


% ADSP Toolbox: Version 2.0 
% For use with "Analog and Digital Signal Processing", 2nd Ed.
% Published by PWS Publishing Co.
%
% Ashok Ambardar, EE Dept. MTU, Houghton, MI 49931, USA
% http://www.ee.mtu/faculty/akambard.html
% e-mail: [e-mail]
% Copyright (c) 1998


if nargin==0,help wbfm,disp('Strike a key to see results of last example')
pause,wbfm('si',10);return,end

ty=ty(1:2);if nargin<3,m=fix(beta);end
w=sqrt(2*beta);
m0=0;
if rem(m,1)~=0,m0=m;m=2*ceil(beta);end
if nargout==0,n1=max([6,2*ceil(beta),m]);else,n1=m;end,
n=0:n1;
if ty=='si'
ti='Sine wave';s=0;for k=0:n1;s(k+1)=besjn(k,beta);end
elseif ty=='st'
ti='Sawtooth';p=n/w+w/4;q=n/w-w/4;
s=exp(-j*2*pi*p.*p).*(cx(2*p)-cx(2*q)+j*sx(2*p)-j*sx(2*q))/w;
elseif ty=='tr'
ti='Triangular';p=n/w+w/2;q=w/2-n/w;
s=(cos(.5*pi*p.*p).*(cx(p)+cx(q))+sin(.5*pi*p.*p).*(sx(p)+sx(q)))/w;
elseif ty=='sq'
ti='Square wave';s=.5*sinc(.5*(beta-n))+.5*((-1).^n).*sinc(.5*(beta+n));
else,error('Ungecognized modulation type'),return
end
s=abs(s);s0=s(1);s1=s(2:m+1);
cs=s0*s0+2*cumsum(s1.*s1);yi=[s0*s0 cs]';
if m0==0;
sm=cs(m);y=sm;
else,
i=find(cs>=m0);m=i(1);y=m;
sm=cs(m);yi=yi(1:m+1);
end
if nargout==0,
ms=int2str(m);
s=[s(n1+1:-1:2) s];n=-n1:n1;
dtplot(n,s,'.')
bt=num2str(beta);ss=num2str(100*sm);
title([ti ' modulation, b = ' bt '  %pwr in ' ms ' harmonics = ' ss '%']),
end