数字信号处理(基于计算机的方法 第三版)matlab程序

% Program 2_1

% Generation of the ensemble average

%

R = 50;

m = 0:R-1;

s = 2*m.*(0.9.^m); % Generate the uncorrupted signal

d = rand(R,1)-0.5; % Generat

e the random noise

x1 = s+d';

stem(m,d);

xlabel('Time index n');ylabel('Amplitude'); title('Noise');

pause

for n = 1:50;

d = rand(R,1)-0.5;

x = s + d';

x1 = x1 + x;

end

x1 = x1/50;

stem(m,x1);

xlabel('Time index n');ylabel('Amplitude'); title('Ensemble average');

% Program 2_2

% Generation of complex exponential sequence

%

a = input('Type in real exponent = ');

b = input('Type in imaginary exponent = ');

c = a + b*i;

K = input('Type in the gain constant = ');

N = input ('Type in length of sequence = ');

n = 1:N;

x = K*exp(c*n);%Generate the sequence

stem(n,real(x));%Plot the real part

xlabel('Time index n');ylabel('Amplitude');

title('Real part');

disp('PRESS RETURN for imaginary part');

pause

stem(n,imag(x));%Plot the imaginary part

xlabel('Time index n');ylabel('Amplitude');

title('Imaginary part');

% Program 2_3

% Generation of real exponential sequence

%

a = input('Type in argument = ');

K = input('Type in the gain constant = ');

N = input ('Type in length of sequence = '); n = 0:N;

x = K*a.^n;

stem(n,x);

xlabel('Time index n');ylabel('Amplitude'); title(['\alpha = ',num2str(a)]);

% Program 2_4

% Signal Smoothing by a Moving-Average Filter %

R = 50;

d = rand(R,1)-0.5;

m = 0:1:R-1;

s = 2*m.*(0.9.^m);

x = s + d';

plot(m,d,'r-',m,s,'b--',m,x,'g:')

xlabel('Time index n'); ylabel('Amplitude') legend('d[n]','s[n]','x[n]');

pause

M = input('Number of input samples = ');

b = ones(M,1)/M;

y = filter(b,1,x);

plot(m,s,'r-',m,y,'b--')

legend('s[n]','y[n]');

xlabel ('Time index n');ylabel('Amplitude')

% Program 2_5

% Illustration of Median Filtering

%

N = input('Median Filter Length = ');

R = 50; a = rand(1,R)-0.4;

b = round(a); % Generate impulse noise

m = 0:R-1;

s = 2*m.*(0.9.^m); % Generate signal

x = s + b; % Impulse noise corrupted signal y = medfilt1(x,3); % Median filtering

subplot(2,1,1)

stem(m,x);axis([0 50 -1 8]);

xlabel('n');ylabel('Amplitude');

title('Impulse Noise Corrupted Signal');

subplot(2,1,2)

stem(m,y);

xlabel('n');ylabel('Amplitude');

title('Output of Median Filter');

% Program 2_6

% Illustration of Convolution

%

a = input('Type in the first sequence = ');

b = input('Type in the second sequence = ');

c = conv(a, b);

M = length(c)-1;

n = 0:1:M;

disp('output sequence =');disp(c)

stem(n,c)

xlabel('Time index n'); ylabel('Amplitude');

% Program 2_7

% Computation of Cross-correlation Sequence

%

x = input('Type in the reference sequence = ');

y = input('Type in the second sequence = ');

% Compute the correlation sequence

n1 = length(y)-1; n2 = length(x)-1;

r = conv(x,fliplr(y));

k = (-n1):n2';

stem(k,r);

xlabel('Lag index'); ylabel('Amplitude');

v = axis;

axis([-n1 n2 v(3:end)]);

% Program 2_8

% Computation of Autocorrelation of a

% Noise Corrupted Sinusoidal Sequence

%

N = 96;

n = 1:N;

x = cos(pi*0.25*n); % Generate the sinusoidal sequence

d = rand(1,N) - 0.5; % Generat

e the noise sequence

y = x + d; % Generate the noise-corrupted sinusoidal sequence r = conv(y, fliplr(y)); % Compute the correlation sequence

k = -28:28;

stem(k, r(68:124));

xlabel('Lag index'); ylabel('Amplitude');

% Program 3_1

% Discrete-Time Fourier Transform Computation

%

% Read in the desired number of frequency samples

k = input('Number of frequency points = ');

% Read in the numerator and denominator coefficients num = input('Numerator coefficients = ');

den = input('Denominator coefficients = ');

% Compute the frequency response

w = 0:pi/(k-1):pi;

h = freqz(num, den, w);

% Plot the frequency response

subplot(2,2,1)

plot(w/pi,real(h));grid

title('Real part')

xlabel('\omega/\pi'); ylabel('Amplitude')

subplot(2,2,2)

plot(w/pi,imag(h));grid

title('Imaginary part')

xlabel('\omega/\pi'); ylabel('Amplitude')

subplot(2,2,3)

plot(w/pi,abs(h));grid

title('Magnitude Spectrum')

xlabel('\omega/\pi'); ylabel('Magnitude')

subplot(2,2,4)

plot(w/pi,angle(h));grid

title('Phase Spectrum')

xlabel('\omega/\pi'); ylabel('Phase, radians')

% Program 3_2

% Generate the filter coefficients

h1 = ones(1,5)/5; h2 = ones(1,14)/14;

% Compute the frequency responses

[H1,w] = freqz(h1, 1, 256);

[H2,w] = freqz(h2, 1, 256);

% Compute and plot the magnitude responses

m1 = abs(H1); m2 = abs(H2);

plot(w/pi,m1,'r-',w/pi,m2,'b--');

ylabel('Magnitude'); xlabel('\omega/\pi');

legend('M=5','M=14');

pause

% Compute and plot the phase responses

ph1 = angle(H1)*180/pi; ph2 = angle(H2)*180/pi;

plot(w/pi,ph1,w/pi,ph2);

ylabel('Phase, degrees');xlabel('\omega/\pi');

legend('M=5','M=14');

% Program 3_3

% Set up the filter coefficients

b = [-6.76195 13.456335 -6.76195];

% Set initial conditions to zero values

zi = [0 0];

% Generate the two sinusoidal sequences

n = 0:99;

x1 = cos(0.1*n);

x2 = cos(0.4*n);

% Generate the filter output sequence

y = filter(b, 1, x1+x2, zi);

% Plot the input and the output sequences

plot(n,y,'r-',n,x2,'b--',n,x1,'g-.');grid

axis([0 100 -1.2 4]);

ylabel('Amplitude'); xlabel('Time index n');

legend('y[n]','x2[n]','x1[n]')

% Program 4_1

% 4-th Order Analog Butterworth Lowpass Filter Design %

format long

% Determine zeros and poles

[z,p,k] = buttap(4);

disp('Poles are at');disp(p);

% Determine transfer function coefficients

[pz, pp] = zp2tf(z, p, k);

% Print coefficients in descending powers of s

disp('Numerator polynomial'); disp(pz)

disp('Denominator polynomial'); disp(real(pp))

omega = [0: 0.01: 5];

% Compute and plot frequency response

h = freqs(pz,pp,omega);

plot (omega,20*log10(abs(h)));grid

xlabel('Normalized frequency'); ylabel('Gain, dB');

% Program 4_2

% Program to Design Butterworth Lowpass Filter

%

% Type in the filter order and passband edge frequency N = input('Type in filter order = ');

Wn = input('3-dB cutoff angular frequency = ');

% Determine the transfer function

[num,den] = butter(N,Wn,'s');

% Compute and plot the frequency response

omega = [0: 200: 12000*pi];

h = freqs(num,den,omega);

plot (omega/(2*pi),20*log10(abs(h)));

xlabel('Frequency, Hz'); ylabel('Gain, dB');

% Program 4_3

% Program to Design Type 1 Chebyshev Lowpass Filter

%

% Read in the filter order, passband edge frequency

% and passband ripple in dB

N = input('Order = ');

Fp = input('Passband edge frequency in Hz = ');

Rp = input('Passband ripple in dB = ');

% Determine the coefficients of the transfer function [num,den] = cheby1(N,Rp,2*pi*Fp,'s');

% Compute and plot the frequency response

omega = [0: 200: 12000*pi];

h = freqs(num,den,omega);

plot (omega/(2*pi),20*log10(abs(h)));

xlabel('Frequency, Hz'); ylabel('Gain, dB');

% Program 4_4

% Program to Design Elliptic Lowpass Filter

%

% Read in the filter order, passband edge frequency, % passband ripple in dB and minimum stopband

% attenuation in dB

N = input('Order = ');

Fp = input('Passband edge frequency in Hz = ');

Rp = input('Passband ripple in dB = ');

Rs = input('Minimum stopband attenuation in dB = '); % Determine the coefficients of the transfer function [num,den] = ellip(N,Rp,Rs,2*pi*Fp,'s');

% Compute and plot the frequency response

omega = [0: 200: 12000*pi];

h = freqs(num,den,omega);

plot (omega/(2*pi),20*log10(abs(h)));

xlabel('Frequency, Hz'); ylabel('Gain, dB');

function y = circonv(x1, x2)

% Develops a sequence y obtained by the circular

% convolution of two equal-length sequences x1 and x2 L1 = length(x1);

L2 = length(x2);

if L1 ~= L2, error('Sequences of unequal length'), end y = zeros(1,L1);

x2tr = [x2(1) x2(L2:-1:2)];

for k = 1:L1,

sh = circshift(x2tr', k-1)';

h = x1.*sh;

disp(sh);

y(k) = sum(h);

end

function [Y] = haar_1D(X)

% Computes the Haar transform of the input vector X

% The length of X must be a power-of-2

% Recursively builds the Haar matrix H

v = log2(length(X)) - 1;

H = [1 1;1 -1];

for m = 1:v,

A = [kron(H,[1 1]);

2^(m/2).*kron(eye(2^m),[1 -1])];

H = A;

end

Y = H*double(X(:));

% Program 5_1

% Illustration of DFT Computation

%

% Read in the length N of sequence and the desired

% length M of the DFT

N = input('Type in the length of the sequence = '); M = input('Type in the length of the DFT = ');

% Generate the length-N time-domain sequence

u = [ones(1,N)];

% Compute its M-point DFT

U = fft(u,M);

% Plot the time-domain sequence and its DFT

t = 0:1:N-1;

stem(t,u)

title('Original time-domain sequence')

xlabel('Time index n'); ylabel('Amplitude')

pause

subplot(2,1,1)

k = 0:1:M-1;

stem(k,abs(U))

title('Magnitude of the DFT samples')

xlabel('Frequency index k'); ylabel('Magnitude') subplot(2,1,2)

stem(k,angle(U))

title('Phase of the DFT samples')

xlabel('Frequency index k'); ylabel('Phase')

% Program 5_2

% Illustration of IDFT Computation

%

% Read in the length K of the DFT and the desired % length N of the IDFT

K = input('Type in the length of the DFT = ');

N = input('Type in the length of the IDFT = ');

% Generate the length-K DFT sequence

k = 0:K-1;

V = k/K;

% Compute its N-point IDFT

v = ifft(V,N);

% Plot the DFT and its IDFT

stem(k,V)

xlabel('Frequency index k'); ylabel('Amplitude') title('Original DFT samples')

pause

subplot(2,1,1)

n = 0:N-1;

stem(n,real(v))

title('Real part of the time-domain samples')

xlabel('Time index n'); ylabel('Amplitude')

subplot(2,1,2)

stem(n,imag(v))

title('Imaginary part of the time-domain samples')

xlabel('Time index n'); ylabel('Amplitude')

% Program 5_3

% Numerical Computation of Fourier transform Using DFT

k = 0:15; w = 0:511;

x = cos(2*pi*k*3/16);% Generate the length-16 sinusoidal sequence X = fft(x); % Compute its 16-point DFT

XE = fft(x,512); % Compute its 512-point DFT

% Plot the frequency response and the 16-point DFT samples

plot(k/16,abs(X),'o', w/512,abs(XE))

xlabel('\omega/\pi'); ylabel('Magnitude')

% Program 5_4

% Linear Convolution Via the DFT

%

% Read in the two sequences

x = input('Type in the first sequence = ');

h = input('Type in the second sequence = ');

% Determine the length of the result of convolution

L = length(x)+length(h)-1;

% Compute the DFTs by zero-padding

XE = fft(x,L); HE = fft(h,L);

% Determine the IDFT of the product

y1 = ifft(XE.*HE);

% Plot the sequence generated by DFT-based convolution and

% the error from direct linear convolution

n = 0:L-1;

subplot(2,1,1)

stem(n,y1)

xlabel('Time index n');ylabel('Amplitude');

title('Result of DFT-based linear convolution')

y2 = conv(x,h);

error = y1-y2;

subplot(2,1,2)

stem(n,error)

xlabel('Time index n');ylabel('Amplitude')

title('Error sequence')

% Program 5_5

% Illustration of Overlap-Add Method

%

% Generate the noise sequence

R = 64;

d = rand(R,1)-0.5;

% Generate the uncorrupted sequence and add noise

k = 0:R-1;

s = 2*k.*(0.9.^k);

x = s + d';

% Read in the length of the moving average filter

M = input('Length of moving average filter = ');

% Generate the moving average filter coefficients

h = ones(1,M)/M;

% Perform the overlap-add filtering operation

y = fftfilt(h,x,4);

% Plot the results

plot(k,s,'r-',k,y,'b--')

xlabel('Time index n');ylabel('Amplitude')

legend('s[n]','y[n]')

function Factors = factorize(polyn)

format long; Factors = [];

% Use threshold of 1e-8 instead of 0 to account for

% precision effects

THRESH = 1e-8;

%

proots = roots(polyn); % get the zeroes of the polynomial

len = length(proots); % get the number of zeroes

%

while(len > 1)

if(abs(imag(proots(1))) < THRESH) % if the zero is a real zero fac = [1 -real(proots(1))];

% construct the factor with proots(1) as zero

Factors = [Factors;[fac 0]];

else % if the zero has imaginary part get all zeroes whose

% imag part is -ve of imaginary part of proots(1)

negimag = imag(proots)+imag(proots(1));

% get all zeroes which have same real part as proot(1)

samereal = real(proots)-real(proots(1));

%find the complex conjugate zero

index = find(abs(negimag)

fac = [1 -2*real(proots(1)) abs(proots(1))^2];

%form 2nd order factor

Factors = [Factors;fac];

else % if the complex conjugate does not exist

fac = [1 -proots(1)];

Factors = [Factors;[fac 0]];

end

end

polyn = deconv(polyn,fac);

%deconvolve the 1st/2nd order factor from polyn

proots = roots(polyn); %determine the new zeros

len = length(polyn); %determine the number of zeros

end

% Program 6_1

% Determination of the Factored Form

% of a Rational z-Transform

%

num = input('Type in the numerator coefficients = ');

den = input('Type in the denominator coefficients = ');

K = num(1)/den(1);

Numfactors = factorize(num)

Denfactors = factorize(den)

disp('Numerator factors');disp(Numfactors);

disp('Denominator factors');disp(Denfactors);

disp('Gain constant');disp(K);

zplane(num,den)

% Program 6_2

% Determination of the Rational z-Transform

% from its Poles and Zeros

%

format long

zr = input('Type in the zeros as a row vector = ');

pr = input('Type in the poles as a row vector = ');

% Transpose zero and pole row vectors

z = zr'; p = pr';

k = input('Type in the gain constant = ');

[num, den] = zp2tf(z, p, k);

disp('Numerator polynomial coefficients'); disp(num);

disp('Denominator polynomial coefficients'); disp(den);

%Program 6_3

% Partial-Fraction Expansion of Rational z-Transform

%

num = input('Type in numerator coefficients = ');

den = input('Type in denominator coefficients = ');

[r,p,k] = residuez(num,den);

disp('Residues');disp(r')

disp('Poles');disp(p')

disp('Constants');disp(k)

% Program 6_4

% Partial-Fraction Expansion to Rational z-Transform

%

r = input('Type in the residues = ');

p = input('Type in the poles = ');

k = input('Type in the constants = ');

[num, den] = residuez(r,p,k);

disp('Numerator polynomial coefficients'); disp(num)

disp('Denominator polynomial coefficients'); disp(den)

% Program 6_5

% Power Series Expansion of a Rational z-Transform

%

% Read in the number of inverse z-transform coefficients to be computed L = input('Type in the length of output vector = ');

% Read in the numerator and denominator coefficients of

% the z-transform

num = input('Type in the numerator coefficients = ');

den = input('Type in the denominator coefficients = ');

% Compute the desired number of inverse transform coefficients

[y,t] = impz(num,den,L);

disp('Coefficients of the power series expansion');

disp(y')

% Program 6_6

% Power Series Expansion of a Rational z-Transform

%

% Read in the number of inverse z-transform coefficients

% to be computed

N = input('Type in the length of output vector = ');

% Read in the numerator and denominator coefficients of % the z-transform

num = input('Type in the numerator coefficients = '); den = input('Type in the denominator coefficients = '); % Compute the desired number of inverse transform

% coefficients

x = [1 zeros(1, N-1)];

y = filter(num, den, x);

disp('Coefficients of the power series expansion'); disp(y)

%Program_6_7

% Stability Test

%

num = input('Type in the numerator vector = ');

den = input('Type in the denominator vector = ');

N = max(length(num),length(den));

x = 1; y0 = 0; S = 0;zi = zeros(1,N-1);

for n = 1:1000

[y,zf] = filter(num,den,x,zi);

if abs(y) < 0.000001, break, end

x = 0;

S = S + abs(y);

y0 = y;zi = zf;

end

if n < 1000

disp('Stable Transfer Function');

else

disp('Unstable Transfer Function');

end

% Program 7_1

% Illustration of Deconvolution

%

Y = input('Type in the convolved sequence = ');

H = input('Type in the convolving sequence = ');

[X,R] = deconv(Y,H);

disp('Sequence x[n]');disp(X);

disp('Remainder Sequence r[n]');disp(R);

% Program 7_2

% Stability Test of a Rational Transfer Function

%

% Read in the polynomial coefficients

den = input('Type in the denominator coefficients ='); % Generate the stability test parameters

k = poly2rc(den);

knew = fliplr(k');

disp('The stability test parameters are');disp(knew); stable = all(abs(k) < 1)

% Program 8_2

% Factorization of a Rational IIR Transfer Function

%

format short

num = input('Numerator coefficients = ');

den = input('Denominator coefficients = '); Numfactors = factorize(num);

Denfactors = factorize(den);

disp('Numerator Factors'),disp(Numfactors)

disp('Denominator Factors'),disp(Denfactors)

% Program 8_3

% Parallel Realizations of an IIR Transfer Function

%

num = input('Numerator coefficient vector = ');

den = input('Denominator coefficient vector = ');

[r1,p1,k1] = residuez(num,den);

[r2,p2,k2] = residue(num,den);

disp('Parallel Form I')

disp('Residues are');disp(r1);

disp('Poles are at');disp(p1);

disp('Constant value');disp(k1);

disp('Parallel Form II')

disp('Residues are');disp(r2);

disp('Poles are at');disp(p2);

disp('Constant value');disp(k2);

% Program 8_4

% Cascaded Lattice Realization of an

% Allpass Transfer Function

%

format long

den = input('Denominator coefficient vector = ');

k = poly2rc(den);

knew = fliplr(k);

disp('The lattice section multipliers are');disp(knew');

% Program 8_5

% Realization of Gray-Markel Cascaded Lattice Structure %

% den is the denominator coefficient vector

% num is the numerator coefficient vector

% k is the lattice parameter vector

% alpha is the vector of feedforward multipliers

%

format long

% Read in the transfer function coefficients

num = input('Numerator coefficient vector = ');

den = input('Denominator coefficient vector = ');

num = num/den(1);

den = den/den(1);

[k,alpha] = tf2latc(num,den);

disp('Lattice parameters are');disp(k');

disp('Feedforward multipliers are');disp(fliplr(alpha'));

% Program 8_6

% Transfer Function of Gray-Markel Cascaded

% Lattice Structure from the Lattice and

% Feedforward Parameters

% k1 is the lattice parameter vector

% alpha is the vector of feedforward multipliers

% den is the denominator coefficient vector

% num is the numerator coefficient vector

%

format long

% Read in the lattice and feedforward parameters

k1 = input('Lattice parameter vector = ');

alpha = input('Feedforward parameter vector = '); [num,den] = latc2tf(k1,fliplr(alpha));

disp('Numerator coefficients are');disp(num)

disp('Denominator coefficients are');disp(den)

% Program 9_1

% Elliptic IIR Lowpass Filter Design

%

Wp = input('Normalized passband edge = ');

Ws = input('Normalized stopband edge = ');

Rp = input('Passband ripple in dB = ');

Rs = input('Minimum stopband attenuation in dB = '); [N,Wn] = ellipord(Wp,Ws,Rp,Rs)

[b,a] = ellip(N,Rp,Rs,Wn);

[h,omega] = freqz(b,a,256);

plot (omega/pi,20*log10(abs(h)));grid;

xlabel('\omega/\pi'); ylabel('Gain, dB');

title('IIR Elliptic Lowpass Filter');

% Program 9_2

% Type 1 Chebyshev IIR Highpass Filter Design

%

Wp = input('Normalized passband edge = ');

Ws = input('Normalized stopband edge = ');

Rp = input('Passband ripple in dB = ');

Rs = input('Minimum stopband attenuation in dB = '); [N,Wn] = cheb1ord(Wp,Ws,Rp,Rs);

[b,a] = cheby1(N,Rp,Wn,'high');

[h,omega] = freqz(b,a,256);

plot (omega/pi,20*log10(abs(h)));grid;

xlabel('\omega/\pi'); ylabel('Gain, dB');

title('Type I Chebyshev Highpass Filter');

% Program 9_3

% Design of IIR Butterworth Bandpass Filter

%

Wp = input('Passband edge frequencies = ');

Ws = input('Stopband edge frequencies = ');

Rp = input('Passband ripple in dB = ');

Rs = input('Minimum stopband attenuation = ');

[N,Wn] = buttord(Wp, Ws, Rp, Rs);

[b,a] = butter(N,Wn);

[h,omega] = freqz(b,a,256);

gain = 20*log10(abs(h));

plot (omega/pi,gain);grid;

xlabel('\omega/\pi'); ylabel('Gain, dB');

title('IIR Butterworth Bandpass Filter');

% Program 9_4

% Group-delay equalization of an IIR filter.

%

[n,d] = ellip(4,1,35,0.3);

[GdH,w] = grpdelay(n,d,512);

plot(w/pi,GdH); grid

xlabel('\omega/\pi'); ylabel('Group delay, samples'); title('Original Filter');

pause

F = 0:0.001:0.3;

g = grpdelay(n,d,F,2); % Equalize the passband

Gd = max(g)-g;

% Design the allpass delay equalizer

[num,den,tau] = iirgrpdelay(8, F, [0 0.3], Gd); [GdA,w] = grpdelay(num,den,512);

plot(w/pi,GdH+GdA); grid

xlabel('\omega/\pi');ylabel('Group delay, samples'); title('Group Delay Equalized Filter');

《应用计算方法教程》matlab作业二

6-1 试验目的计算特征值，实现算法 试验容：随机产生一个10阶整数矩阵，各数均在-5和5之间。 （1） 用MATLAB 函数“eig ”求矩阵全部特征值。 （2） 用幂法求A 的主特征值及对应的特征向量。 （3） 用基本QR 算法求全部特征值（可用MATLAB 函数“qr ”实现矩阵的QR 分解）。 原理 幂法：设矩阵A 的特征值为12n ||>||||λλλ≥???≥并设A 有完全的特征向量系12,,,n χχχ???(它们线性无关)，则对任意一个非零向量0n V R ∈所构造的向量序列1k k V AV -=有11()lim ()k j k k j V V λ→∞ -=， 其中()k j V 表示向量的第j 个分量。 为避免逐次迭代向量k V 不为零的分量变得很大（1||1λ>时）或很小（1||1λ<时），将每一步的k V 按其模最大的元素进行归一化。具体过程如下： 选择初始向量0V ，令1max(),,,1k k k k k k k V m V U V AU k m +===≥，当k 充分大时1111,max()max() k k U V χλχ+≈ ≈。 QR 法求全部特征值： 111 11222 111 ,1,2,3,k k k k k A A Q R R Q A Q R k R Q A Q R +++==????==??=???? ??????==?? 由于此题的矩阵是10阶的，上述算法计算时间过长，考虑采用改进算法——移位加速。迭 代格式如下： 1 k k k k k k k k A q I Q R A R Q q I +-=?? =+? 计算k A 右下角的二阶矩阵() () 1,1 1,() (),1 ,k k n n n n k k n n n n a a a a ----?? ? ??? 的特征值()()1,k k n n λλ-，当()()1,k k n n λλ-为实数时，选k q 为()()1,k k n n λλ-中最接近(),k n n a 的。 程序

数字信号处理Matlab实现实例(推荐给学生)

数字信号处理Matlab 实现实例 第１章离散时间信号与系统 例1-1 用MATLAB计算序列{-2 0 1 –1 3}和序列{1 2 0 -1}的离散卷积。 解 MATLAB程序如下： a=[-2 0 1 -1 3]; b=[1 2 0 -1]; c=conv(a,b); M=length(c)-1; n=0:1:M; stem(n,c); xlabel('n'); ylabel('幅度'); 图1.1给出了卷积结果的图形，求得的结果存放在数组c中为：{-2 -4 1 3 1 5 1 -3}。 例1-2 用MATLAB计算差分方程 当输入序列为时的输出结果。 解 MATLAB程序如下： N=41; a=[0.8 -0.44 0.36 0.22]; b=[1 0.7 -0.45 -0.6]; x=[1 zeros(1,N-1)];

k=0:1:N-1; y=filter(a,b,x); stem(k,y) xlabel('n');ylabel('幅度') 图 1.2 给出了该差分方程的前41个样点的输出，即该系统的单位脉冲响应。 例1-3 用MATLAB 计算例1-2差分方程 所对应的系统函数的DTFT 。 解 例1-2差分方程所对应的系统函数为： 123 123 0.80.440.360.02()10.70.450.6z z z H z z z z -------++= +-- 其DTFT 为 23230.80.440.360.02()10.70.450.6j j j j j j j e e e H e e e e ωωωω ωωω--------++= +-- 用MATLAB 计算的程序如下： k=256; num=[0.8 -0.44 0.36 0.02]; den=[1 0.7 -0.45 -0.6]; w=0:pi/k:pi; h=freqz(num,den,w); subplot(2,2,1); plot(w/pi,real(h));grid title('实部') xlabel('\omega/\pi');ylabel('幅度')

数字信号处理MATLAB中FFT实现

MATLAB中FFT的使用方法 说明：以下资源来源于《数字信号处理的MATLAB实现》万永革主编 一.调用方法 X=FFT(x)； X=FFT(x，N)； x=IFFT(X); x=IFFT(X,N) 用MATLAB进行谱分析时注意： （1）函数FFT返回值的数据结构具有对称性。 例： N=8; n=0:N-1; xn=[43267890]; Xk=fft(xn) → Xk= 39.0000-10.7782+6.2929i0-5.0000i 4.7782-7.7071i 5.0000 4.7782+7.7071i0+5.0000i-10.7782-6.2929i Xk与xn的维数相同，共有8个元素。Xk的第一个数对应于直流分量，即频率值为0。 （2）做FFT分析时，幅值大小与FFT选择的点数有关，但不影响分析结果。在IFFT时已经做了处理。要得到真实的振幅值的大小，只要将得到的变换后结果乘以2除以N即可。 二.FFT应用举例 例1：x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)。采样频率fs=100Hz，分别绘制N=128、1024点幅频图。

clf; fs=100;N=128;%采样频率和数据点数 n=0:N-1;t=n/fs;%时间序列 x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);%信号 y=fft(x,N);%对信号进行快速Fourier变换 mag=abs(y);%求得Fourier变换后的振幅 f=n*fs/N;%频率序列 subplot(2,2,1),plot(f,mag);%绘出随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=128');grid on; subplot(2,2,2),plot(f(1:N/2),mag(1:N/2));%绘出Nyquist频率之前随频率变化的振幅xlabel('频率/Hz'); ylabel('振幅');title('N=128');grid on; %对信号采样数据为1024点的处理 fs=100;N=1024;n=0:N-1;t=n/fs; x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);%信号 y=fft(x,N);%对信号进行快速Fourier变换 mag=abs(y);%求取Fourier变换的振幅 f=n*fs/N; subplot(2,2,3),plot(f,mag);%绘出随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=1024');grid on; subplot(2,2,4) plot(f(1:N/2),mag(1:N/2));%绘出Nyquist频率之前随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=1024');grid on; 运行结果：

计算方法_全主元消去法_matlab程序

%求四阶线性方程组的MA TLAB程序 clear Ab=[0.001 2 1 5 1; 3 - 4 0.1 -2 2; 2 -1 2 0.01 3; 1.1 6 2.3 9 4];%增广矩阵 num=[1 2 3 4];%未知量x的对应序号 for i=1:3 A=abs(Ab(i:4,i:4));%系数矩阵取绝对值 [r,c]=find(A==max(A(:))); r=r+i-1;%最大值对应行号 c=c+i-1;%最大值对应列号 q=Ab(r,:),Ab(r,:)=Ab(i,:),Ab(i,:)=q;%行变换 w=Ab(:,c),Ab(:,c)=Ab(:,i),Ab(:,i)=w;%列变换 n=num(i),num(i)=num(c),num(c)=n;%列变换引起未知量x次序变化for j=i:3 Ab(j+1,:)=-Ab(j+1,i)*Ab(i,:)/Ab(i,i)+Ab(j+1,:);%消去过程 end end %最后得到系数矩阵为上三角矩阵 %回代算法求解上三角形方程组 x(4)=Ab(4,5)/Ab(4,4); x(3)=(Ab(3,5)-Ab(3,4)*x(4))/Ab(3,3); x(2)=(Ab(2,5)-Ab(2,3)*x(3)-Ab(2,4)*x(4))/Ab(2,2); x(1)=(Ab(1,5)-Ab(1,2)*x(2)-Ab(1,3)*x(3)-Ab(1,4)*x(4))/Ab(1,1); for s=1:4 fprintf('未知量x%g =%g\n',num(s),x(s)) end %验证如下 %A=[0.001 2 1 5 1; 3 -4 0.1 -2 2;2 -1 2 0.01 3; 1.1 6 2.3 9 4]; %b=[1 2 3 4]'; %x=A\b; %x1= 1.0308 %x2= 0.3144 %x3= 0.6267 %x4= -0.0513

实验一 基于Matlab的数字信号处理基本

实验一 基于Matlab 的数字信号处理基本操作 一、 实验目的：学会运用MA TLAB 表示的常用离散时间信号；学会运用MA TLAB 实现离 散时间信号的基本运算。 二、 实验仪器：电脑一台，MATLAB6.5或更高级版本软件一套。 三、 实验内容： （一） 离散时间信号在MATLAB 中的表示 离散时间信号是指在离散时刻才有定义的信号，简称离散信号，或者序列。离散序列通常用)(n x 来表示，自变量必须是整数。 离散时间信号的波形绘制在MATLAB 中一般用stem 函数。stem 函数的基本用法和plot 函数一样，它绘制的波形图的每个样本点上有一个小圆圈，默认是空心的。如果要实心，需使用参数“fill ”、“filled ”，或者参数“.”。由于MATLAB 中矩阵元素的个数有限，所以MA TLAB 只能表示一定时间范围内有限长度的序列；而对于无限序列，也只能在一定时间范围内表示出来。类似于连续时间信号，离散时间信号也有一些典型的离散时间信号。 1. 单位取样序列 单位取样序列)(n δ，也称为单位冲激序列，定义为 ) 0() 0(0 1)(≠=?? ?=n n n δ 要注意，单位冲激序列不是单位冲激函数的简单离散抽样，它在n =0处是取确定的值1。在MATLAB 中，冲激序列可以通过编写以下的impDT .m 文件来实现，即 function y=impDT(n) y=(n==0); %当参数为0时冲激为1,否则为0 调用该函数时n 必须为整数或整数向量。 【实例1-1】 利用MATLAB 的impDT 函数绘出单位冲激序列的波形图。 解：MATLAB 源程序为 >>n=-3:3; >>x=impDT(n); >>stem(n,x,'fill'),xlabel('n'),grid on >>title('单位冲激序列') >>axis([-3 3 -0.1 1.1]) 程序运行结果如图1-1所示。 图1-1 单位冲激序列

青岛理工大学临沂年数字信号处理及MATLAB试卷

A卷

一、[15分] 1、10 2、f>=2fh

3、()()()y n x n h n =* 4、1 -az -11a 或者-z z ，a 1 -z 或1-1-az -1z 5、对称性 、 可约性 、 周期性 6、191点，256 7、典范型、级联型、并联型 8、T ω = Ω，)2 tan(2ω T = Ω或)2arctan(2T Ω=ω。 二、[20分] 1、C 2、 A 3、 C 4、C 5、B 6、D 7、B 8、A 9、D 10、A (CACCB DBADA) 三、[15分] 1、（5分） 混叠失真：不满足抽样定理的要求。 改善方法：增加记录长度 频谱泄漏：对时域截短，使频谱变宽拖尾，称为泄漏 改善方法：1）增加w (n )长度 2）缓慢截短 栅栏效应：DFT 只计算离散点（基频F0的整数倍处）的频谱，而不是连续函数。 改善方法：增加频域抽样点数N （时域补零），使谱线更密 2、（5分） 3、 （5分） IIR 滤波器: 1）系统的单位抽样相应h (n )无限长 2）系统函数H (z )在有限z 平面（ ）上有极点存在 3）存在输出到输入的反馈，递归型结构 Fir 滤波器: ? 1）系统的单位冲激响应h (n )在有限个n 处不为零； ? 2）系统函数 在||0 z >处收敛，在 处只有零点，即有限z 平面只有零点，而全部极点都在z =0处； ? 3）机构上主要是非递归结构，没有输入到输出的反馈，但有些结构中也包含有反馈的递归部分。 四、计算题（40分） 1、（12分）解： 解： 对上式两边取Z 变换，得： ()H z ||0z >

数字信号处理指导书matlab版

实验1 时域离散信号的产生 一、实验目的 学会运用MATLAB 产生常用离散时间信号。 二、实验涉及的matlab 子函数 1、square 功能：产生矩形波 调用格式： x=square(t)；类似于sin （t ），产生周期为2*pi ，幅值为+—1的方波。 x=square(t ，duty)；产生制定周期的矩形波，其中duty 用于指定脉冲宽度与整个周期的比例。 2、rand 功能：产生rand 随机信号。 调用格式： x=rand （n ，m ）；用于产生一组具有n 行m 列的随机信号。 三、实验原理 在时间轴的离散点上取值的信号，称为离散时间信号。通常，离散时间信号用x （n ）表示，其幅度可以在某一范围内连续取值。 由于信号处理所用的设备主要是计算机或专用的信号处理芯片，均以有限的位数来表示信号的幅度，因此，信号的幅度也必须“量化”，即取离散值。我们把时间和幅度上均取离散值的信号称为时域离散信号或数字信号。 在MATLAB 中，时域离散信号可以通过编写程序直接生成，也可以通过对连续信号的等间隔抽样获得。 下面介绍常用的时域离散信号及其程序。 1、单位抽样序列 ? ? ?≠==000 1)(k k k δ MATLAB 源程序为

1) function [x,n] = impuls (n0,n1,n2) % Generates x(n) = delta(n-n0); n=n0 处建立一个单位抽样序列% [x,n] = impuls (n0,n1,n2) if ((n0 < n1) | (n0 > n2) | (n1 > n2)) error('arguments must satisfy n1 <= n0 <= n2') end n = [n1:n2]; x = [zeros(1,(n0-n1)), 1, zeros(1,(n2-n0))]; 将上述文件存为：impuls.m，在命令窗口输入 n0=0,n1=-10,n2=11; [x,n]=impuls (n0,n1,n2); stem(n,x,’filled’) 2)n1=-5;n2=5;n0=0; n=n1:n2; x=[n==n0]; stem(n,x,'filled','k'); axis([n1,n2,1.1*min(x),1.1*max(x)]); title('单位脉冲序列'); xlabel('时间(n)'); ylabel('幅度x(n)'); 3)n1=-5;n2=5;k=0; n=n1:n2; nt=length(n); %求n点的个数 nk=abs(k-n1)+1; %确定k在n序列中的位置 x=zeros(1,nt); %对所有样点置0 x(nk)=1; %对抽样点置1 stem(n,x,'filled','k'); axis([n1,n2,0,1.1*max(x)]); title('单位脉冲序列'); xlabel('时间(n)'); Ylabel('幅度x(n)');

王能超 计算方法——算法设计及MATLAB实现课后代码

第一章插值方法 1.1Lagrange插值 1.2逐步插值 1.3分段三次Hermite插值 1.4分段三次样条插值 第二章数值积分 2.1 Simpson公式 2.2 变步长梯形法 2.3 Romberg加速算法 2.4 三点Gauss公式 第三章常微分方程德差分方法 3.1 改进的Euler方法 3.2 四阶Runge-Kutta方法 3.3 二阶Adams预报校正系统 3.4 改进的四阶Adams预报校正系统 第四章方程求根 4.1 二分法 4.2 开方法 4.3 Newton下山法 4.4 快速弦截法 第五章线性方程组的迭代法 5.1 Jacobi迭代 5.2 Gauss-Seidel迭代 5.3 超松弛迭代 5.4 对称超松弛迭代 第六章线性方程组的直接法 6.1 追赶法 6.2 Cholesky方法 6.3 矩阵分解方法 6.4 Gauss列主元消去法

第一章插值方法 1.1Lagrange插值 计算Lagrange插值多项式在x=x0处的值. MATLAB文件：（文件名：Lagrange_eval.m）function [y0,N]= Lagrange_eval(X,Y,x0) %X,Y是已知插值点坐标 %x0是插值点 %y0是Lagrange插值多项式在x0处的值 %N是Lagrange插值函数的权系数 m=length(X); N=zeros(m,1); y0=0; for i=1:m N(i)=1; for j=1:m if j~=i; N(i)=N(i)*(x0-X(j))/(X(i)-X(j)); end end y0=y0+Y(i)*N(i); end 用法》X=[…];Y=[…]; 》x0= ; 》[y0,N]= Lagrange_eval(X,Y,x0) 1.2逐步插值 计算逐步插值多项式在x=x0处的值. MATLAB文件：（文件名：Neville_eval.m）function y0=Neville_eval(X,Y,x0) %X,Y是已知插值点坐标 %x0是插值点 %y0是Neville逐步插值多项式在x0处的值 m=length(X); P=zeros(m,1); P1=zeros(m,1); P=Y; for i=1:m P1=P; k=1; for j=i+1:m k=k+1;

数字信号处理的MATLAB实现

昆明理工大学信息工程与自动化学院学生实验报告 （2011—2012 学年第二学期） 课程名称：数字信号处理开课实验室：信自楼111 2012 年 5 月 31 日年级、专业、班生医学号姓 名 成绩 实验项目名称数字信号处理的matlab 实现指导教师 教 师 评语教师签名： 年月日 一．实验目的 熟练掌握matlab的基本操作。 了解数字信号处理的MATLAB实现。 二．实验设备 安装有matlab的PC机一台。 三．实验内容 .1.求信号x(n)=cos(6.3Пn/3)+cos(9.7Пn/30)+cos(15.3Пn/30),0≤n≤29的幅度频谱. 2. 用冲击响应不变法设计一个Butterworth低通数字滤波器,要求参数为: Wp=0.2Пαp=1dB Ws=0.3Пαs=15dB 3.用双线性变换法设计一个Chebyshev高通IIR滤波器,要求参数为: Wp=0.6Пαp=1dB Ws=0.4586Пαs=15dB 4.用窗函数法设计一个低通FIR滤波器,要求参数为: Wp=0.2Пαp=0.3dB Ws=0.25Пαs=50dB 5.用频率抽样法设计一个带通FIR滤波器,要求参数为: W1s=0.2П W1p=0.35П W2p=0.65П W2s=0.8П αs=60dB αp=1dB 6.根据 4 点矩形序列，( n ) = [1 1 1 1] 。做 DTFT 变换，再做 4 点 DFT 变换。然后分别补零做 8 点 DFT 及 16 点 DFT。 7.调用filter解差分方程，由系统对u(n)的响应判断稳定性 8编制程序求解下列系统的单位冲激响应和阶跃响应。 y[n]+ 0.75y[n -1]+ 0.125y[n -2] = x[n]- x[n -1] 四．实验源程序 1. n=[0:1:29]; x=cos(6.3*pi*n/30)+cos(9.7*pi*n/30)+cos(15.3*pi*n/30);

(整理)matlab16常用计算方法.

常用计算方法 1．超越方程的求解 一超越方程为 x (2ln x – 3) -100 = 0 求超越方程的解。 [算法]方法一：用迭代算法。将方程改为 01002ln()3 x x =- 其中x 0是一个初始值，由此计算终值x 。取最大误差为e = 10-4，当| x - x 0| > e 时，就用x 的值换成x 0的值，重新进行计算；否则| x - x 0| < e 为止。 [程序]P1_1abs.m 如下。 %超越方程的迭代算法 clear %清除变量 x0=30; %初始值 xx=[]; %空向量 while 1 %无限循环 x=100/(2*log(x0)-3); %迭代运算 xx=[xx,x]; %连接结果 if length(xx)>1000,break ,end %如果项数太多则退出循环（暗示发散） if abs(x0-x)<1e-4,break ,end %当精度足够高时退出循环 x0=x; %替换初值 end %结束循环 figure %创建图形窗口 plot(xx,'.-','LineWidth',2,'MarkerSize',12)%画迭代线'.-'表示每个点用.来表示，再用线连接 grid on %加网格 fs=16; %字体大小 title('超越方程的迭代折线','fontsize',fs)%标题 xlabel('\itn','fontsize',fs) %x 标签 ylabel('\itx','fontsize',fs) %y 标签 text(length(xx),xx(end),num2str(xx(end)),'fontsize',fs)%显示结果 [图示]用下标作为自变量画迭代的折线。如P0_20_1图所示，当最大误差为10-4时，需要迭代19次才能达到精度，超越方程的解为27.539。 [算法]方法二：用求零函数和求解函数。将方程改为函数 100()2ln()3f x x x =-- MATLAB 求零函数为fzero ，fzero 函数的格式之一是 x = fzero(f,x0) 其中，f 表示求解的函数文件，x0是估计值。fzero 函数的格式之二是 x = fzero(f,[x1,x2])

数字信号处理MATLAB实验1

实验一熟悉MATLAB环境 一、实验目的 (1)熟悉MATLAB的主要操作命令。 (2)学会简单的矩阵输入和数据读写。 (3)掌握简单的绘图命令。 (4)用MATLAB编程并学会创建函数。 (5)观察离散系统的频率响应。 二、实验内容 认真阅读本章附录，在MATLAB环境下重新做一遍附录中的例子，体会各条命令的含义。在熟悉了MATLAB基本命令的基础上，完成以下实验。 上机实验内容： (1)数组的加、减、乘、除和乘方运算。输入A=[1234]，B=[345 6]，求C=A+B，D=A-B，E=A.*B，F=A./B，G=A.^B并用stem语句画出 A、B、C、D、E、F、G。 (2)用MATLAB实现以下序列。 a）x(n)=0.8n0≤n≤15 b）x(n)=e(0.2+3j)n0≤n≤15 c）x(n)=3cos(0.125πn+0.2π)+2sin(0.25πn+0.1π)0≤n≤15 (n)=x(n+16),绘出四个d）将c)中的x(n)扩展为以16为周期的函数x 16 周期。 (n)=x(n+10),绘出四个e）将c)中的x(n)扩展为以10为周期的函数x 10 周期。

(3)x(n)=[1,-1,3,5]，产生并绘出下列序列的样本。 a）x 1(n)=2x(n+2)-x(n-1)-2x(n) b）∑=-=5 1k 2) k n (nx (n) x (4)绘出下列时间函数的图形，对x轴、y轴以及图形上方均须加上适当的标注。 a)x(t)=sin(2πt)0≤t≤10s b)x(t)=cos(100πt)sin(πt) 0≤t≤4s (5)编写函数stepshift(n0,n1,n2)实现u(n-n0),n1

数字信号处理基本知识点Matlab实现

数字信号处理（第二版） 绪论 1.4 MATLAB 在信号处理中的应用简介 MATLAB 是美国Mathworks 公司于1984年推出的一套高性能的数值计算和可视化软件，它集数值分析、矩阵运算、信号处理、系统仿真和图形显示于一体，从而被广泛地应用于科学计算、控制系统、信息处理等领域的分析、仿真和设计工作。 MATLAB 软件包括五大通用功能：数值计算功能（Numeric ），符号运算功能（Symbolic ）；数据可视化功能（Graphic ），数据图形文字统一处理功能（Notebook ）和建模仿真可视化功能（Simulink ）。该软件有三大特点：一是功能强大；二是界面友善、语言自然；三是开放性强。目前，Mathworks 公司已推出30多个应用工具箱。MA TLAB 在线性代数、矩阵分析、数值及优化、数理统计和随机信号分析、电路与系统、系统动力学、信号和图像处理、控制理论分析和系统设计、过程控制、建模和仿真、通信系统、以及财政金融等众多领域的理论研究和工程设计中得到了广泛应用。 2.10 离散时间信号与系统的Matlab 表示 2.10.1 离散时间信号的表示和运算 1、基本序列的Matlab 表示 单位采样序列 在MA TLAB 中，单位采样序列可以通过编写以下的DTimpulse .m 文件来实现，即 function y=DTimpulse (n) y=(n==0); %当参数为0时冲激为1,否则为0 调用该函数时n 必须为整数或整数向量。 单位阶跃序列 在MA TLAB 中，单位阶跃序列可以通过编写DTu .m 文件来实现，即 function y=DTu (n) y=n>=0; %当参数为非负时输出1 调用该函数时n 必须为整数或整数向量。 矩形序列 用MA TLAB 表示矩形序列可根据公式()()()N R n u n u n N =--并利用DTu 函数生成，即 function y=DTR(n,N) y=DTu(n)-DTu(n-N); 调用该函数时n 必须为整数或整数向量，N 必须为整数。 实指数序列 用MA TLAB 表示实指数序列()(),n x n a u n n N a R =∈∈，即

南京理工大学数字信号处理matlab上机完美版

1.已知3阶椭圆IIR数字低通滤波器的性能指标为：通带截止频率0.4π，通带波纹为0.6dB，最小阻带衰减为32dB。设计一个6阶全通滤波器对其通带的群延时进行均衡。绘制低通滤波器和级联滤波器的群延时。 %Q1_solution %ellip(N,Ap,Ast,Wp) %N--->The order of the filter %Ap-->ripple in the passband %Ast->a stopband Rs dB down from the peak value in the passband %Wp-->the passband width [be,ae]=ellip(3,0.6,32,0.4); hellip=dfilt.df2(be,ae); f=0:0.001:0.4; g=grpdelay(hellip,f,2); g1=max(g)-g; [b,a,tau]=iirgrpdelay(6,f,[0 0.4],g1); hallpass=dfilt.df2(b,a); hoverall=cascade(hallpass,hellip); hFVT=fvtool([hellip,hoverall]); set(hFVT,'Filter',[hellip,hoverall]); legend(hFVT,'Lowpass Elliptic filter','Compensated filter'); clear; [num1,den1]=ellip(3,0.6,32,0.4); [GdH,w]=grpdelay(num1,den1,512); plot(w/pi,GdH); grid xlabel('\omega/\pi'); ylabel('Group delay, samples'); F=0:0.001:0.4; g=grpdelay(num1,den1,F,2); % Equalize the passband Gd=max(g)-g; % Design the allpass delay equalizer [num2,den2]=iirgrpdelay(6,F,[0,0.4],Gd); [GdA,w] = grpdelay(num2,den2,512); hold on; plot(w/pi,GdH+GdA,'r');

matlab用于计算方法的源程序

1、Newdon迭代法求解非线性方程 function [x k t]=NewdonToEquation(f,df,x0,eps) %牛顿迭代法解线性方程 %[x k t]=NewdonToEquation(f,df,x0,eps) %x:近似解 %k:迭代次数 %t:运算时间 %f:原函数，定义为内联函数 ?:函数的倒数，定义为内联函数 %x0:初始值 %eps:误差限 % %应用举例： %f=inline('x^3+4*x^2-10'); ?=inline('3*x^2+8*x'); %x=NewdonToEquation(f,df,1,0.5e-6) %[x k]=NewdonToEquation(f,df,1,0.5e-6) %[x k t]=NewdonToEquation(f,df,1,0.5e-6) %函数的最后一个参数也可以不写。默认情况下，eps=0.5e-6 %[x k t]=NewdonToEquation(f,df,1) if nargin==3 eps="0".5e-6; end tic; k=0; while 1 x="x0-f"(x0)./df(x0); k="k"+1; if abs(x-x0) < eps || k >30 break; end x0=x; end t=toc; if k >= 30 disp('迭代次数太多。'); x="0"; t="0"; end

2、Newdon迭代法求解非线性方程组 function y="NewdonF"(x) %牛顿迭代法解非线性方程组的测试函数 %定义是必须定义为列向量 y(1,1)=x(1).^2-10*x(1)+x(2).^2+8; y(2,1)=x(1).*x(2).^2+x(1)-10*x(2)+8; return; function y="NewdonDF"(x) %牛顿迭代法解非线性方程组的测试函数的导数 y(1,1)=2*x(1)-10; y(1,2)=2*x(2); y(2,1)=x(2).^+1; y(2,2)=2*x(1).*x(2)-10; return; 以上两个函数仅供下面程序的测试 function [x k t]=NewdonToEquations(f,df,x0,eps) %牛顿迭代法解非线性方程组 %[x k t]=NewdonToEquations(f,df,x0,eps) %x:近似解 %k:迭代次数 %t:运算时间 %f:方程组（事先定义） ?:方程组的导数（事先定义） %x0:初始值 %eps:误差限 % %说明：由于虚参f和df的类型都是函数，使用前需要事先在当前目录下采用函数M文件定义% 另外在使用此函数求解非线性方程组时，需要在函数名前加符号“@”，如下所示 % %应用举例： %x0=[0,0];eps=0.5e-6; %x=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps) %[x k]=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps) %[x k t]=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps) %函数的最后一个参数也可以不写。默认情况下，eps=0.5e-6 %[x k t]=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps)

基于MATLAB的数字信号处理

数字信号处理课程设计报告题目：语音数字信号处理与分析及 Matlab实现 系别通信工程 专业班级 学生姓名 学号 指导教师 提交日期

摘要 本次课程设计综合利用数字信号处理的理论知识进行语音信号的频谱分析，通过理论推导得出相应结论，再利用MATLAB作为编程工具进行计算机实现，从而加深对所学知识的理解，建立概念。本次课程设计要求利用MATLAB对语音信号进行分析和处理，要求学生采集语音信号后，在MATLAB软件平台进行频谱分析；并对所采集的语音信号加入干扰噪声，对加入噪声的信号进行频谱分析，设计合适的滤波器滤除噪声，恢复原信号。待处理语音信号是一个在20Hz～20kHz 频段的低频信号。采用了高效快捷的开发工具——MATLAB，实现了语音信号的采集，对语音信号加噪声及设计滤波器滤除噪声的一系列工作。利用采样原理设计了高通滤波器、低通滤波器、带通滤波器、带阻滤波器。同学通过查阅资料自己获得程序进行滤波器的设计，能过得到很好的锻炼。 关键词：MATLAB滤波器数字信号处理

目录 第一章绪论 (1) 1.1设计的目的及意义 (1) 1.2设计要求 (1) 1.3设计内容 (1) 第二章系统方案论证 (3) 2.1设计方案分析 (3) 2.2实验原理 (3) 第三章信号频谱分析 (6) 3.1原始信号及频谱分析 (6) 3.2加入干扰噪声后的信号及频谱分析 (7) 第四章数字滤波器的设计与实现 (11) 4.1高通滤波器的设计 (11) 4.2低通滤波器的设计 (12) 4.3带通滤波器的设计 (15) 4.4带阻滤波器的设计 (16) 第五章课程设计总结 (19) 参考文献 (20) 附录Ⅰ..................................................................................I 附录Ⅱ................................................................................II

0计算方法及MATLAB实现简明讲义课件PPS8-1欧拉龙格法

第8章 常微分方程初值问题数值解法 8.1 引言 8.2 欧拉方法 8.3 龙格-库塔方法 8.4 单步法的收敛性与稳定性 8.5 线性多步法

8.1 引 言 考虑一阶常微分方程的初值问题 00(,),[,],(). y f x y x a b y x y '=∈=（1.1） （1.2） 如果存在实数 ，使得 121212(,)(,).,R f x y f x y L y y y y -≤-?∈（1.3） 则称 关于 满足李普希茨(Lipschitz ）条件， 称为 的李普希茨常数（简称Lips.常数）. 0>L f y L f (参阅教材386页)

计算方法及MATLAB 实现 所谓数值解法，就是寻求解 在一系列离散节点 )(x y <<<<<+121n n x x x x 上的近似值 . ,,,,,121+n n y y y y 相邻两个节点的间距 称为步长. n n n x x h -=+1 如不特别说明，总是假定 为定数， ),2,1( ==i h h i 这时节点为 . ) ,2,1,0(0 =+=i nh x x n 初值问题(1.1),(1.2)的数值解法的基本特点是采取 “步进式”. 即求解过程顺着节点排列的次序一步一步地向前推进. 00(,),[,], (). y f x y x a b y x y '=∈=

描述这类算法，只要给出用已知信息 ,,,21--n n n y y y 计算 的递推公式. 1+n y 一类是计算 时只用到前一点的值 ，称为单步法. 1+n y n y 另一类是用到 前面 点的值 ， 1+n y k 11,,,+--k n n n y y y 称为 步法. k 其次，要研究公式的局部截断误差和阶，数值解 与 精确解 的误差估计及收敛性，还有递推公式的计算 稳定性等问题. n y )(n x y 首先对方程 离散化，建立求数值解的递推 公式. ),(y x f y ='

数字信号处理实验 matlab版 线性相位FIR数字滤波器

实验23 线性相位FIR数字滤波器 (完美格式版，本人自己完成，所有语句正确，不排除极个别错误，特别适用于山大，勿用冰点等工具下载，否则下载之后的word格式会让很多部分格式错误，谢谢) XXXX学号姓名处XXXX 一、实验目的 1 加深对线性相位FIR数字滤波器特性的理解。 2 掌握线性相位滤波器符幅特性和零极点分布的研究方法。 3 了解用MATLAB研究线性相位滤波器特性时程序编写的思路和方法。 二、实验内容 1 线性相位FIR滤波器的特性 2 第一类线性相位滤波器(类型Ⅰ) 3 第二类线性相位滤波器(类型Ⅱ) 4 第三类线性相位滤波器(类型Ⅲ) 5 第四类线性相位滤波器(类型Ⅳ) 6 线性相位FIR数字滤波器零点分布特点 三、实验环境 MATLAB7.0 四、实验原理 1.线性相位FIR滤波器的特性 与IIR滤波器相比，FIR滤波器在保证幅度特性满足技术要求的同时，很容易做到有严格的线性相位特性。设FIR滤波器单位脉冲响应h(n)长度为N，其系统函数为 ∑-=- = 1 N n n z h(n) H(z) 当滤波器的系数N满足一定的对称条件时，就可以获得线性相位。线性相位FIR滤波器共分为四种类型，分别为： (1)类型Ⅰ，系数对称，即h(n)＝h(N－1－n)，N为奇数。 (2)类型Ⅱ，系数对称，即h(n)＝h(N－1－n)，N为偶数。 (3)类型Ⅲ，系数反对称，即h(n)＝－h(N－1－n)，N为奇数。 (4)类型Ⅳ，系数反对称，即h(n)＝－h(N－1－n)，N为偶数。 对于上述四类线性相位FIR滤波器，参考文献［1］中提供了一段通用程序，对考虑正负号的幅度频率特性(简称符幅特性)进行求解，程序名为amplres.m，程序如下：function[A,w,type,tao]=amplres(h) N=length(h);tao=(N-1)/2; L=floor((N-1)/2); n=1:L+1; w=[0:500]*2*pi/500; if all(abs(h(n)-h(N-n+1))<1e-10)

数字信号处理MATLAB仿真(DOC)

实验一 数字信号处理的Matlab 仿真 一、实验目的 1、掌握连续信号及其MATLAB 实现方法； 2、掌握离散信号及其MATLAB 实现方法 3、掌握离散信号的基本运算方法，以及MATLAB 实现 4、了解离散傅里叶变换的MATLAB 实现 5、了解IIR 数字滤波器设计 6、了解FIR 数字滤波器设计1 二、实验设备 计算机，Matlab 软件 三、实验内容 （一）、 连续信号及其MATLAB 实现 1、 单位冲击信号 ()0,0()1,0t t t dt εεδδε-?=≠??=?>??? 例1.1：t=1/A=50时，单位脉冲序列的MATLAB 实现程序如下： clear all; t1=-0.5:0.001:0; A=50; A1=1/A; n1=length(t1); u1=zeros(1,n1); t2=0:0.001:A1;

t0=0; u2=A*stepfun(t2,t0); t3=A1:0.001:1; n3=length(t3); u3=zeros(1,n3); t=[t1 t2 t3]; u=[u1 u2 u3]; plot(t,u) axis([-0.5 1 0 A+2]) 2、 任意函数 ()()()f t f t d τδττ+∞ -∞=-? 例1.2：用MATLAB 画出如下表达式的脉冲序列 ()0.4(2)0.8(1) 1.2() 1.5(1) 1.0(2)0.7(3)f n n n n n n n δδδδδδ=-+-+++++++ clear all; t=-2:1:3; N=length(t); x=zeros(1,N); x(1)=0.4; x(2)=0.8 x(3)=1.2; x(4)=1.5; x(5)=1.0;

用MATLAB实现结构可靠度计算.

用MATLAB实现结构可靠度计算 口徐华…朝泽刚‘u刘勇‘21 。 (【l】中国地质大学(武汉工程学院湖北?武汉430074; 12】河海大学土木工程学院江苏?南京210098 摘要:Matlab提供了各种矩阵的运算和操作,其中包含结构可靠度计算中常用的各种数值计算方法工具箱,本文从基本原理和相关算例分析两方面,阐述利用Matlab,编制了计算结构可靠度Matlab程.序,使得Matlab-语言在可靠度计算中得到应用。 关键词:结构可靠度Matlab软件最优化法 中图分类号:TP39文献标识码:A文章编号:1007-3973(200902-095-Ol 1结构可靠度的计算方法 当川概率描述结构的可靠性时,计算结构可靠度就是计算结构在规定时问内、规定条件F结构能够完成预定功能的概率。 从简单到复杂或精确稃度的不同,先后提出的可靠度计算方法有一次二阶矩方法、二次二阶矩方法、蒙特卡洛方法以及其他方法。一次■阶矩方法又分为。I-心点法和验算点法,其中验算点法足H前可靠度分析最常川的方法。 2最优化方法计算可靠度指标数学模型 由结构111n个任意分布的独立随机变量一,x:…以表示的结构极限状态方程为:Z=g(■.托…t=0,采用R-F将非正念变量当罱正态化,得到等效正态分布的均值o:和标准差虹及可靠度指标B,由可靠度指标B的几何意义知。o;辟

开始时验算点未知,把6看成极限状态曲面上点P(■,爿:---37,的函数,通过优化求解,找到B最小值。求解可靠皮指标aJ以归结为以下约束优化模型: rain睁喜t华,2 s.,.Z=g(工i,x2’,…,工:=0 如极限状态方栉巾某个变最(X。可用其他变量表示,则上述模型jfIJ‘转化为无约束优化模型: 。。B!:手f生丛r+阻:坚:坠:盐尘}二剐 t∞oY?’【叫,J 3用MATLAB实现结构可靠度计算 3.1Matlab简介 Matlab是++种功能强、效率高、便.丁.进行科学和工程计算的交互式软件包,汇集了人量数学、统计、科学和工程所需的函数,MATI.AB具有编程简甲直观、用户界mf友善、开放性强等特点。将MATLAB用于蒙特卡罗法的一个显著优点是它拥有功能强大的随机数发生器指令。 3.2算例 3.2.I例:已知非线形极限状态方程z=g(t r'H=567f r-0.5H2=0’f、r服从正态分布。IIf=0.6,o r=0.0786;la|_ 2.18,o r_0.0654;H服从对数正态分布。u H= 3218,O。 =0.984。f、r、H相互独立,求可靠度指标B及验算点(,,r’,H‘。 解:先将H当量正念化:h=ln H服从正态分布,且 ,‘-““了:等专虿’=,。49?口二-、『五ir面_。。3