matlab数字信号处理源码

简单画一维函数曲线

clear all;
     close all;
     t= 0:0.01:10*pi;
     A=1;
     a=-0.4;
     ft=cos(a*t);
     st=sin(a*t);
     axis([-40,5,-2,2]);
     plot(t,ft);
      title('f(t)=cos(a*t)');
      
subplot(2,1,1);
plot(t,ft);
xlabel('w');
ylabel('F1(jw)');
title(' ft=cos(a*t)');

subplot(2,1,2);
plot(t,st);
xlabel('w');
ylabel('F1(jw)');
title(' st=sin(a*t)');

画幅频曲线

b0 = 0.05634;
b1 = [1  1];
b2 = [1 -1.0166 1];
a1 = [1 -0.683];
a2 = [1 -1.4461 0.7957];

b = b0*conv(b1,b2);
a = conv(a1,a2);

[h,w] = freqz(b,a,'whole',2001);


plot(w/pi,20*log10(abs(h)))
ax = gca;
ax.YLim = [-100 20];
ax.XTick = 0:.5:2;
xlabel('Normalized Frequency (\times\pi rad/sample)')
ylabel('Magnitude (dB)')

画信号和频谱

离散信号

Fs = 1000;            % Sampling frequency
T = 1/Fs;             % Sampling period
L = 10;             % Length of signal
t = (0:L-1)*T;        % Time vector
S =  sin(2*pi*120*t);

plot(t,S)
title('Signal')
xlabel('t (milliseconds)')
ylabel('X(t)')

DFS

%{使用单个输入参数 x 的 fft 计算输入向量或矩阵的 DFT。如果 x 是向量，fft 计算向量的 DFT；如果 x 是矩形数组，fft 计算每个数组列的 DFT。
%例如，创建时间向量和信号：
%}
t = 0:1/100:1-1/100;                     % Time vector
x = sin(2*pi*5*t) + sin(2*pi*40*t);      % Signal
%计算信号的 DFT 以及变换后的序列的幅值和相位。通过将小幅值变换值设置为零来减少计算相位时的舍入误差。

y = fft(x);                               % Compute DFT of x
m = abs(y);                               % Magnitude
y(m<1e-6) = 0;
p = unwrap(angle(y));                     % Phase
%要以度为单位绘制幅值和相位，请键入以下命令：

f = (0:length(y)-1)*100/length(y);        % Frequency vector%归为0~100的量

subplot(3,1,1)
plot(t,x)
xlabel('n');
ylabel('x(n)');
title('Signal')


subplot(3,1,2)
plot(f,m)
xlabel('w');
ylabel('|F1(jw)|');
title('Magnitude')
ax = gca;
%ax.XTick = [15 40 60 85];

subplot(3,1,3)
plot(f,p*180/pi)
xlabel('w');
ylabel('arg(F1(jw))');
title('Phase')
%ax = gca;
%ax.XTick = [15 40 60 85];

FFT

F=input('输入信号频率f/Hz：');                   %输入信号频率f%
N=input('输入采样点数N：');                      %输入采样点数N%
T=input('输入采样周期T/s：');                    %输入采样间隔T%
%采样
for j=0:1:N-1
x(j+1)=sin(2*pi*f*j*T);               
end

%补0
insert=input('是否补0？[0-否；1-是]:');    
if insert==1                              
    ZERO=input('补零数：');    
for j=N:1:N+ZERO-1
x(j+1)=0;                      
end
    N=N+ZERO;                           
end

%码位倒置%
M=log2(N);              
for t=1:1:N                           
    s=dec2bin(t-1,M);                 
    s=fliplr(s);                   
    s=bin2dec(s);                 
    A(s+1)=x(t);             
end

%按时间抽取的FFT蝶形运算%
for L=1:1:M                 
for J=0:1:(2^(L-1)-1)          
for k=(J+1):2^L:N       
           T=A(k)+A(k+2^(L-1))*exp((-i*2*pi*J*2^(M-L))/N);
A(k+2^(L-1))=A(k)-A(k+2^(L-1))*exp((-i*2*pi*J*2^(M-L))/N);
A(k)=T;             
end
end
end

%模值归一化%
x=abs(A);                           
y=max(x);                            
X=x/y;

%绘图
for j=1:1:N
stem(j-1,X(j));                 
hold on
end
axis([0 N 0 1]);