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\');
请发表评论