MATLAB的一些应用--最近用的比较多
1、MATLAB分析信号的频谱
快速Fourier变换(FFT)是离散傅里叶变换的快速算法,他是根据离散傅里叶变换的奇、偶、虚、实等特性,对离散傅里叶变换的算法进行改进获得的。
针对几个个简单例子介绍一下：
（1、）假设数据采集频率为1000Hz,一个信号包含频率为50Hz、振幅为0.7的正弦波和频率为120Hz、振幅为1的正弦波,噪声为零平均值的随机噪声
用FFT方法分析其频谱方法Matlab程序如下:
clear Fs = 1000; % 采样频率 T = 1/Fs; % 采样时间 L = 1000; % 信号长度 t = (0:L-1)*T; % 时间向量 x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t); y = x + 2*randn(size(t)); % 加噪声正弦信号 figure(1) plot(Fs*t(1:50),y(1:50)) title('零平均值噪音信号'); xlabel('time (milliseconds)') NFFT = 2^nextpow2(L); % Next power of 2 from length of y Y = fft(y,NFFT)/L; f = Fs/2*linspace(0,1,NFFT/2); figure(2) plot(f,2*abs(Y(1:NFFT/2))) title('y(t)单边振幅频谱') xlabel('Frequency (Hz)') ylabel('|Y(f)|')
结果如下：
（2）产生余弦信号以作频谱分析：余弦信号y=cos(2π*f*t)；信号频率为f=10Hz；时宽:1s 采样率为fs=100Hz；
MATLAB程序：
clear all; f=10; fs=100; T=1; n=round(T*fs);%采样点个数 t=linspace(0,T,n); y=cos(2*pi*f/fs*[0:n-1]); figure(1) plot(t,y); title('余弦信号时域'); xlabel('t/s'); ylabel('幅度'); %用fft函数对产生的余弦信号作频谱分析： %注意：该步骤得到的是0~fs内的频谱。 fft_y=fft(y); f=linspace(0,fs,n); figure(2) plot(f,abs(fft_y)); title('余弦信号频谱(fft)'); xlabel('f/Hz'); ylabel('幅度'); % 用fftshift函数得到-fs/2~fs/2内的频谱： fftshift_y=fftshift(fft_y); f=linspace(-fs/2,fs/2,n); figure(3) plot(f,abs(fftshift_y)); title('余弦信号频谱FFTshift'); xlabel('f/Hz'); ylabel('幅度');
结果：
可以看到10Hz处有峰值，90Hz的峰值是-10Hz的峰值向右频谱搬移fs=100Hz得到的。
由于实信号频谱幅度关于原点对称，可以看到10Hz与-10Hz处的两个峰值。
2、MATLAB中几种采样方法及实现
X(t)的时域信号
syms x;
>> f=sym('cos(2/3*pi*x)');
>> ezplot(f,[0,40]);
采样信号及频域波形
w=-pi:0.01*pi:pi; n=0:40; x=cos(2/3*pi*n); X=x*exp(-j*n'*w); subplot(211); stem(n,x,'filled'); xlabel('n'); title('x[n]'); subplot(212); plot(w/pi,X);
过采样：
w=-pi:0.01*pi:pi; n=0:40; x=cos(2/3*pi*n); X=x*exp(-j*n'*w); subplot(311); stem(n,x,'filled'); xlabel('n'); title('x[n]'); subplot(312); plot(w/pi,abs(X)); xlabel('Omega/pi'); title('Magnitude of X'); subplot(313); plot(w/pi,angle(X)); xlabel('Omega/pi'); title('Phase of X');
临界采样：
w=-pi:0.01*pi:pi; n=0:1.5:40; x=cos(2/3*pi*n); X=x*exp(-j*n'*w); subplot(311); stem(n,x,'filled'); xlabel('n'); title('x[n]'); subplot(312); plot(w/pi,abs(X)); xlabel('Omega/pi'); title('Magnitude of X'); subplot(313); plot(w/pi,angle(X)); xlabel('Omega/pi'); title('Phase of X');
欠采样：
w=-pi:0.01*pi:pi; n=0:3:40; x=cos(2/3*pi*n); X=x*exp(-j*n'*w); subplot(311); stem(n,x,'filled'); xlabel('n'); title('x[n]'); subplot(312); plot(w/pi,abs(X)); xlabel('Omega/pi'); title('Magnitude of X'); subplot(313); plot(w/pi,angle(X)); xlabel('Omega/pi'); title('Phase of X');
信号重构--临界采样
n=-100:100; t=-30:0.005:30; wc=2; Ts=pi/wc; ws=2*pi/Ts; m=n*Ts; f=sinc(m/pi); ft=f*Ts*wc*sinc((wc/pi)*(ones(length(m),1)*t-m'*ones(1,length(t))))./pi; t1=-9*pi:Ts:9*pi; f1=sinc(t1/pi); subplot(311); stem(t1,f1,'filled'); xlabel('kTs'); ylabel('kTs'); title('临界采样信号'); subplot(312); plot(t,ft); title('临界采样信号重构信号'); xlabel('t'); ylabel('f(t)'); subplot(313); plot(t,ft-sinc(t/pi)); title('重构信号与原信号误差'); xlabel('t');
信号重构--欠采样
n=-100:100; t=-30:0.005:30; wc=0.5; Ts=pi/wc; ws=2*pi/Ts; m=n*Ts; f=sinc(m/pi); ft=f*Ts*wc*sinc((wc/pi)*(ones(length(m),1)*t-m'*ones(1,length(t))))./pi; t1=-9*pi:Ts:9*pi; f1=sinc(t1/pi); subplot(311); stem(t1,f1,'filled'); xlabel('kTs'); ylabel('kTs'); title('欠采样信号'); subplot(312); plot(t,ft); xlabel('t'); ylabel('f(t)'); title('欠采样信号重构信号'); subplot(313); plot(t,ft-sinc(t/pi)); title('重构信号与原信号误差'); xlabel('t');
信号重构--过采样
n=-100:100; t=-30:0.005:30; wc=4; Ts=pi/wc; ws=2*pi/Ts; m=n*Ts; f=sinc(m/pi); ft=f*Ts*wc*sinc((wc/pi)*(ones(length(m),1)*t-m'*ones(1,length(t))))./pi; t1=-9*pi:Ts:9*pi; f1=sinc(t1/pi); subplot(311); stem(t1,f1,'filled'); xlabel('kTs'); ylabel('kTs'); title('过采样信号'); subplot(312); plot(t,ft); xlabel('t'); ylabel('f(t)'); title('过采样信号重构信号'); subplot(313); plot(t,ft-sinc(t/pi)); title('重构信号与原信号误差'); xlabel('t');
