Matlab 信号处理【问答笔记-1】

1. 关于matlab frezq fir1函数

问：输入信号x为频率0.1 Hz和0.3 Hz的等幅正弦波之和，利用fir1函数设计滤波器h，去除0.3 Hz的正弦信号，得到输出信号y。用freqz函数观察滤波器h的频率响应。显示输入信号x和输出信号y的波形。用fft函数对信号做快速傅里叶变换，并显示输入输出信号的幅频曲线。

答：具体步骤如下，精度自行调整；由于0.1Hz与0.3Hz频率差距太小，容易产生频谱泄漏：

clc,clear,close all;
t = -100:1:100;
L = length(t);
fs = 2;
x1 = sin(2*pi*0.1.*t/fs);
x2 = sin(2*pi*0.3.*t/fs);
x3 = x1 + x2;
%时域信号
figure(1)
subplot(311)
plot(t,x1,"LineWidth",1.5)
grid on 
subplot(312)
plot(t,x2,"LineWidth",1.5)
grid on
subplot(313)
plot(t,x3,"LineWidth",1.5)
grid on
%滤波器设计
wn = 0.2;  %截止频率wn为0.2pi(0.2Hz),wn = 2pi*f/fs
N = 60;  %阶数选择
hn = fir1(N-1,wn,boxcar(N));  %10阶FIR低通滤波器
figure(2)
freqz(hn,1);
figure(3)
y = fftfilt(hn,x3);  %经过FIR滤波器后得到的信号
plot(t,y,"LineWidth",1.5)
grid on
%频谱分析
Y = fft(y);  %输出信号的fft
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = fs*(0:(L/2))/L;
plot(f,P1,"LineWidth",1.5) 
title('Single-Sided Amplitude Spectrum of X(t)')
xlabel('f (Hz)')
ylabel('|P1(f)|')

2. 用MATLAB完成傅里叶变换并绘制出变换后的波形

问：
用MATLAB求出傅里叶变换并绘制图像
用什么函数工具，怎么得出图像呀？

答：连续非周期的符号函数sign(t)求傅里叶变换的结果如下：

图像：

在Matlab里绘制的图像实际上都是转成离散时间信号处理后的结果（计算机里没有真正的模拟量），对连续时间信号（模拟信号）采样逼近得到的，这里得出的频谱（DFT得到的频域离散谱是对连续谱的N点等间隔采样）取连续包络就是sign(t)经连续时间傅里叶变换FT后的图像，就不再补充从FS、FT到DFS、DTFT、DFT（FFT）的过程和关系了，具体推导自行看书研究，写在这里会写到吐血的。

clc,clear,close all;
Fs = 1000;                   % Sampling frequency                    
T = 1/Fs;                    % Sampling period       
L = 1500;                    % Length of signal
t = (-1/2*L:1/2*L-1)*T;      % Time vector
 
%时域
y = sign(t);
figure(1)
plot(t,y)
ylim([-2 2])
grid on
title('Time-domain')
 
%fft近似表示ft
Y = fft(y);
P2 = abs(Y/L);
P0 = fftshift(P2);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
figure(2)
subplot(211)
plot(P0) 
title('Frequency-domain')
xlabel('w')
ylabel('|P2(f)|')
%基于 P2 和信号长度 L 计算单侧频谱 P1
subplot(212)
plot(P1) 
title('Frequency-domain')
xlabel('w')
ylabel('|P1(f)|')

3. 该怎么设计butter worth带阻滤波器的各系数

问：

噪音是7000hz，非噪音是0-5000hz，只看左半边，我该怎么设计呢？参见如下设计：


答：
你想要通过带阻滤波器进行滤波，那噪声信号的频率必须在阻带范围内；有用信号为0~5000Hz，因此通带截止频率fp1可以设置为5500Hz，阻带起始频率fs1必小于7000Hz，可以设为6000Hz，阻带截止频率fs2必须大于7000Hz，可以设为8000Hz，通带起始频率fp2可设置为8500Hz，为了保证奈奎斯特条件，采样频率Fs可以设为44100Hz.

clc,clear,close all;
Fs = 44100;
fp1 = 5500;fp2 = 8500;
fs1 = 6000;fs2 = 8000;
wp1 = fp1/Fs*2*pi;
wp2 = fp2/Fs*2*pi;
wp = [wp1,wp2];
ws1 = fs1/Fs*2*pi;
ws2 = fs2/Fs*2*pi;
ws = [ws1,ws2];
Rp = 1;As = 30;
[n,wc] = buttord(wp/pi,ws/pi,Rp,As);
[b,a] = butter(n,wc,'stop');
freqz(b,a);

0.32π*44100/2π约等于7000Hz，可见噪声频率落在阻带范围内。