网站首页轮播图片素材,女装网站建设规划书怎么写,ui设计周末培训学校,西安专业做网站的公司有哪些目录 一、ADC采样的噪声简介1.1 常见的ADC噪声来源 二、信号的时域到频域转换2.1 傅里叶变换巧记傅里叶变换 三、傅里叶变换和滤波算法工程实现3.1 使用Matlab计算信号时域到频域的变换3.2 使用Matlab去除特定频点噪声寻找峰值算噪声频率构建陷波滤波器滤除噪声频点陷波滤波器与… 目录 一、ADC采样的噪声简介1.1 常见的ADC噪声来源 二、信号的时域到频域转换2.1 傅里叶变换巧记傅里叶变换 三、傅里叶变换和滤波算法工程实现3.1 使用Matlab计算信号时域到频域的变换3.2 使用Matlab去除特定频点噪声寻找峰值算噪声频率构建陷波滤波器滤除噪声频点陷波滤波器与滑动平均滤波器对比 3.3 C 实现ADC采样信号滤波 一、ADC采样的噪声简介
经常用ADC采样的朋友都知道采样的数据不是纯净的信号会伴随着干扰。 一般把周期性的干扰叫噪声既然是周期性的一定能算出周期时长T取倒数就是频率f单位Hz。
1.1 常见的ADC噪声来源
电源是噪声的一大来源如果用的是交流电交流电通过整流电路中的二极管元件过滤掉交流电的负半周然后通过滤波电容滤波让交流电趋于平缓形成直流电。交流信号肯定不能滤干净的这就是电源纹波产生的一大原因电源纹波会导致ADC采样信号出现噪声。
二、信号的时域到频域转换
刚接触信号的时候搞不清楚什么是时域什么是频域大学里的《信号与系统》《控制工程与原理》《信号处理》等书上提过但因为没有结合实际使用学过都忘记了。 时域可以理解成横轴是时间t纵轴是其他有幅值变换的物理量(如电压声响亮度温度振幅等)如下图左图所示 频域是横轴为频率大小纵轴为频率分量的强度(或者叫功率密度) 而傅里叶变换就是时域和频域的联结。
2.1 傅里叶变换
从时域变换到频域用的是傅里叶变换傅里叶觉得计算波形相加和分解太麻烦了发明了一种简单的方法让波形相加能像“112”一样直观他的方法就是把波形从时域转换到频域再相加或者分解计算好后再转回时域波形。我们把傅里叶发明的这套方法叫做“傅里叶变换”和“傅里叶逆变换”。 公式如下如果看不懂一是回去复习微积分二是用我下面的方法巧记 F ( ω ) F { f ( t ) } 1 / 2 π ∫ − ∞ ∞ f ( t ) e − i ω t d t ( 傅里叶变换 ) F(\omega) \mathcal{F}\{f(t)\} {1/2\pi} \int_{-\infty}^{\infty} f(t) e^{-i \omega t} \, dt (傅里叶变换) F(ω)F{f(t)}1/2π∫−∞∞f(t)e−iωtdt(傅里叶变换) f ( t ) F − 1 { F ( ω ) } ∫ − ∞ ∞ F ( ω ) e i ω t d ω ( 傅里叶逆变换 ) f(t) \mathcal{F}^{-1}\{F(\omega)\} \int_{-\infty}^{\infty} F(\omega) e^{i \omega t} \, d\omega (傅里叶逆变换) f(t)F−1{F(ω)}∫−∞∞F(ω)eiωtdω(傅里叶逆变换) 也可以用 ω 2 π f \omega 2\pi f ω2πf带入 F ( f ) F { f ( t ) } ∫ − ∞ ∞ f ( t ) e − 2 π i f t d t F(f) \mathcal{F}\{f(t)\} \int_{-\infty}^{\infty} f(t) e^{-2\pi i f t} \, dt F(f)F{f(t)}∫−∞∞f(t)e−2πiftdt f ( t ) F − 1 { F ( f ) } ∫ − ∞ ∞ F ( f ) e 2 π i f t d f f(t) \mathcal{F}^{-1}\{F(f)\} \int_{-\infty}^{\infty} F(f) e^{2\pi i f t} \, df f(t)F−1{F(f)}∫−∞∞F(f)e2πiftdf
巧记傅里叶变换
先看这个动画 傅里叶变换动画 首先看懂公式最里面的这个 e − i ω t e^{-i \omega t} e−iωt这其实就是三角函数在复数域的延申。 已知欧拉公式 e − i ω t c o s ( ω t ) − i s i n ( ω t ) e^{-i \omega t}cos(\omega t)-isin(\omega t) e−iωtcos(ωt)−isin(ωt) e − i ω t e^{-i \omega t} e−iωt可以看成一个距离原地1单位的球在绕着原点以 ω \omega ω角速度旋转从x轴(实数)看过去就是 c o s ( ω t ) cos(\omega t) cos(ωt)从y轴(虚数)看过去就是 s i n ( ω t ) sin(\omega t) sin(ωt)。
然后再看 f ( t ) e − i ω t f(t) e^{-i \omega t} f(t)e−iωt下面动图最上面的蓝线是时域函数f(t)红线就是 e − i ω t e^{-i \omega t} e−iωt的实数轴部分这里有两个变量 ω \omega ω和t把大家看懵了我们细细捋一下下面的横轴就是tt从(0,∞)的变换都已经显示在数轴上了可以看成静态量导致红线变化的是另一个变量 ω \omega ω它才是动态变量。 f ( t ) e − i ω t f(t) e^{-i \omega t} f(t)e−iωt也就是把 f ( t ) f(t) f(t)和 e − i ω t e^{-i \omega t} e−iωt相乘也就是中间图。可以看到当 e − i ω t e^{-i \omega t} e−iωt变化的和 f ( t ) f(t) f(t)波动频率一致的时候乘积结果就是 s i n 2 ( t ) ( 1 − c o s ( 2 t ) ) / 2 sin^2(t) {(1-cos(2t))}/2 sin2(t)(1−cos(2t))/2。
最后经过积分 ∫ − ∞ ∞ f ( t ) e − i ω t d t \int_{-\infty}^{\infty} f(t) e^{-i \omega t} \, dt ∫−∞∞f(t)e−iωtdt得到的下图。注意看这是对时间t进行积分之前说过时间是静态量所以图中横轴已经不是时间t了而是角频率 ω \omega ω公式 F ( ω ) 1 / 2 π ∫ − ∞ ∞ f ( t ) e − i ω t d t F(\omega) {1/2\pi} \int_{-\infty}^{\infty} f(t) e^{-i \omega t} \, dt F(ω)1/2π∫−∞∞f(t)e−iωtdt表示的是积分结果随着频率变化而变化。
最后一张图和时间已经没关系了横轴已经是频率了这样就完成了时域到频域的转化。
PS:为啥要费劲心思把周期信号从时域转频域关键就是周期周期意味着经过一定时间后会重复所以时间反而不重要了重要的是频率。下图的最高点对应的就是原始信号的频率。
(动图参考来源https://blog.csdn.net/YuHWEI/article/details/136678900) 上面的图能解释傅里叶变换但不是真正的傅里叶变换因为我们只计算了实数轴真正的傅里叶变换后正弦信号是在8.3Hz上的一道冲击信号。 同理多个周期性正弦波叠加也可以从频域上看出频率分量。 在频域上做周期性信号的合成与分解是不是就简单多了呢。 通过不同频率的正弦波叠加想要什么波形都能叠加出来。 反过来任何一种波形都可以分解成不同频率的正弦波叠加。
三、傅里叶变换和滤波算法工程实现
3.1 使用Matlab计算信号时域到频域的变换
matlab使用快速傅里叶变换函数(fft)注意fft后的复数域的量还需要abs()一下才能显示到实数轴上。
data xlsread(origin_signal_data.xlsx); % 读取原始波形数据
voltage data(:, 1); N length(voltage); % 数据长度t 2.5; %总采样时间
T t/N; %采样周期总采样时间/总采样次数
Fs 1/T; %采样率是1s内的采样次数
time 0:5:2500-5; %设置横坐标数组
plot(time, voltage, r, LineWidth, 2); % 绘制单边频谱的幅值谱
xlabel(采样时间 (ms)); % 频率轴标签需要知道采样率才能正确标注
ylabel(电压 (mV)); % 幅值轴标签
title(原始时域波形); % 图表标题
grid on;% data_no_dc timeSeries - signal_mean;
Y fft(voltage); % 计算FFT
P2 abs(Y/N); % 双边频谱的幅值并归一化
P1 P2(1:N/21); % 单边频谱的幅值正频率部分
P1(2:end-1) 2*P1(2:end-1); % 因为FFT结果是对称的所以正频率部分的幅值要乘以2除了第一个和最后一个点f (0:N-1)*(1/N); % 频率向量归一化频率
f_actual f*(Fs); % 实际频率需要知道采样率
f_used f_actual(1:N/21); %使用到的只有一半plot(f_used, P1); % 绘制单边频谱的幅值谱
xlabel(频率 (Hz)); % 频率轴标签需要知道采样率才能正确标注
ylabel(幅值);
title(傅里叶变换后频域波形); % 图表标题
grid on; % 显示网格下面左图的原始数据中还有周期性噪声干扰通过傅里叶变换后看到噪声主要集中在60Hz。
注意我这里是去除了0Hz的数据0Hz是信号中的直流分量信号都是直流分量和交流分量的叠加如果傅里叶变换后0Hz有很高的幅值那就是直流分量我们要处理的是噪声肯定是交流分量所以作图时不要把0Hz算进去。 3.2 使用Matlab去除特定频点噪声
寻找峰值算噪声频率
ignore_below_frequency 1; % 1Hz以下不管因为是信号的直流分量% 寻找傅里叶变换后最大幅值所在的频点
validIndices f_used ignore_below_frequency;
filteredP1 P1(validIndices);
filteredF f_used(validIndices); [peakAmplitude, peakIndex] max(filteredP1);
peakFrequency filteredF(peakIndex); % 获得噪声频率构建陷波滤波器滤除噪声频点
构建陷波滤波器设置滤波频点和带宽
% 设置60Hz陷波滤波器
wo peakFrequency / (Fs / 2);
bw wo / 20; % 设置滤波器带宽
notchFilter designfilt(bandstopiir, ...FilterOrder, 2, ...HalfPowerFrequency1, wo - bw/2, ...HalfPowerFrequency2, wo bw/2, ...DesignMethod, butter);% 使用滤波器
filtered_voltage filtfilt(notchFilter, voltage); 陷波滤波器与滑动平均滤波器对比
滑动平均滤波器就是设置一个窗口窗口内包含多个采样频点比如前7个数据为一个窗口取窗口内数据的平均值为当前点的值然后窗口向后滑动滑动平均算法的核心思想是在一个正弦周期内所有点之和为零。 如果选的窗口刚好是噪声周期那就滤除干扰了。 但滑动平均在较低频率下好用高频后窗口大小无法继续缩小不如陷波滤波。
3.3 C 实现ADC采样信号滤波
matlab验证方法工程实现还是要用C这里我要把算法最终应用到嵌入式设备上去所以使用QT的C来部署验证。
QT C需要使用快速傅里叶变换库和矩阵库支持 在.pro文件里加入fftw3库。
LIBS -lfftw3同时下载Eigen库放到项目文件中并加入头文件路径
INCLUDEPATH ./eigen-3.4.0陷波滤波器没有现场的库需要手动写
QVectordouble applyNotchFilter(const QVectordouble voltageData, double sampleRate, double notchFreq, double qualityFactor)
{QVectordouble filteredData(voltageData.size());double w0 2.0 * M_PI * notchFreq / sampleRate;double alpha sin(w0) / (2.0 * qualityFactor);double b0 1.0;double b1 -2.0 * cos(w0);double b2 1.0;double a0 1.0 alpha;double a1 -2.0 * cos(w0);double a2 1.0 - alpha;b0 / a0;b1 / a0;b2 / a0;a1 / a0;a2 / a0;double x1 0.0, x2 0.0; double y1 0.0, y2 0.0; for (int i 0; i voltageData.size(); i) {double x0 voltageData[i];double y0 b0 * x0 b1 * x1 b2 * x2 - a1 * y1 - a2 * y2;filteredData[i] y0;x2 x1;x1 x0;y2 y1;y1 y0;}return filteredData;
}为了像matlab一样能直观的看到滤波效果用QChart作图
void displayFilterResult(const QVectordouble timeData, const QVectordouble voltageData, const QVectordouble filteredData, const double maxFrequency)
{// 创建QLineSeries对象并设置它们的名称QLineSeries *seriesVoltage new QLineSeries();seriesVoltage-setName(原始电压数据); // 设置名称为“原始电压数据”QLineSeries *seriesFiltered new QLineSeries();seriesFiltered-setName(滤波后电压数据); // 设置名称为“滤波后电压数据”// 填充数据到QLineSeries对象for (int i 0; i timeData.size(); i) {seriesVoltage-append(timeData[i], voltageData[i]);seriesFiltered-append(timeData[i], filteredData[i]);}// 创建QChart对象并添加系列与之前相同QChart *chart new QChart();chart-addSeries(seriesVoltage);chart-addSeries(seriesFiltered);// 默认情况下图例是可见的但你可以通过以下方式访问和配置它QLegend *legend chart-legend();legend-setVisible(true); // 通常情况下这是默认的但你可以显式设置legend-setAlignment(Qt::AlignBottom); // 设置图例的位置这里是在底部// 为图表设置标题与之前相同QString frequency QString::number(maxFrequency, f, 1) Hz;chart-setTitle(滤波前后对比——干扰频率 frequency);// 创建QValueAxis对象并设置轴的范围QValueAxis *axisX new QValueAxis();axisX-setLabelFormat(%g s);chart-addAxis(axisX, Qt::AlignBottom);seriesVoltage-attachAxis(axisX);seriesFiltered-attachAxis(axisX);QValueAxis *axisY new QValueAxis();axisY-setLabelFormat(%g V);chart-addAxis(axisY, Qt::AlignLeft);seriesVoltage-attachAxis(axisY);seriesFiltered-attachAxis(axisY);// 为不同的系列设置不同的颜色seriesVoltage-setPen(QPen(QColor(255, 0, 0))); // 红色seriesFiltered-setPen(QPen(QColor(0, 0, 255))); // 蓝色// 创建QChartView对象并设置QChartQChartView *chartView new QChartView(chart);chartView-setRenderHint(QPainter::Antialiasing);// 创建主窗口并显示图表QMainWindow* window new QMainWindow();window-setCentralWidget(chartView);window-resize(1920, 1080);window-show();
}主函数如下
#include QFile
#include QString
#include QVector
#include QDebug
#include QtCharts/QChartView
#include QtCharts/QLineSeries
#include QtCharts/QValueAxis
#include QMainWindow
#include QApplication#include fftw3.h
#include Eigen/Denseusing namespace Eigen;
using namespace std;QT_CHARTS_USE_NAMESPACEint main(int argc, char *argv[]) {QApplication a(argc, argv);QFile inputFile(./source_signal_data.csv);if (!inputFile.open(QIODevice::ReadOnly | QIODevice::Text)) {qDebug() 文件打开失败;return -1;}QVectordouble voltageData;while (!inputFile.atEnd()) {QByteArray line inputFile.readLine();voltageData.push_back(line.toDouble());}inputFile.close();// 配置周期和频率int N voltageData.size();double t 2.5; // 总时间double T t / N;double Fs 1 / T; // 采样频率// 生成时间数据从0s开始每隔0.005s到2.5sQVectordouble timeData;for (double t 0; t (2.5 - 0.005); t 0.005) {timeData.append(t);}// 将信号数据放到实部中VectorXd voltageEigen(N);for (int i 0; i N; i) {voltageEigen[i] voltageData[i];}fftw_complex* in (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * N);fftw_complex* out (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * N);fftw_plan plan fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);for (int i 0; i N; i) {in[i][0] voltageEigen[i]; // 实部in[i][1] 0.0; // 虚部}fftw_execute(plan); //应用傅里叶变换VectorXcd fftResult(N);for (int i 0; i N; i) {fftResult[i] std::complexdouble(out[i][0], out[i][1]);}// 清楚缓存fftw_destroy_plan(plan);fftw_free(in);fftw_free(out);// 制作频谱VectorXd absFFT fftResult.cwiseAbs() / N;VectorXd P1 absFFT.head(N / 2 1);P1.segment(1, P1.size() - 2) * 2;VectorXd f_actual VectorXd::LinSpaced(N / 2 1, 0, Fs / 2);// 滤除直流分量QVectordouble filteredP1;QVectordouble filteredF;for (int i 0; i f_actual.size(); i) {if (f_actual(i) 1) {filteredP1.push_back(P1(i));filteredF.push_back(f_actual(i));}}// 导出到csvQFile outputFile(../FFT_test/filtered_data.csv);if (outputFile.open(QIODevice::WriteOnly | QIODevice::Text)) {QTextStream out(outputFile);out 频率,幅值\n;for (int i 0; i filteredF.size(); i) {out filteredF[i] , filteredP1[i] \n;}outputFile.close();} else {qDebug() 导出失败 ;}// 寻找噪声最大的幅值和对应的频率auto maxIt max_element(filteredP1.begin(), filteredP1.end());double maxAmplitude 0.0;double maxFrequency 0.0;if (maxIt ! filteredP1.end()){size_t maxIndex distance(filteredP1.begin(), maxIt);maxAmplitude *maxIt;maxFrequency filteredF[maxIndex];}else{qDebug() 未找到噪声频率;}double qualityFactor 1; // 设置陷波滤波带宽QVectordouble filteredData applyNotchFilter(voltageData, Fs, maxFrequency, qualityFactor);// 导出到csvQString notchFileName ./notch_filtered_ QString::number(maxFrequency, f, 1) Hz.csv;QFile notchOutputFile(notchFileName);if (notchOutputFile.open(QIODevice::WriteOnly | QIODevice::Text)) {QTextStream out(notchOutputFile);for (int i 0; i filteredData.size(); i) {out filteredData[i] \n;}notchOutputFile.close();} else {qDebug() 导出失败;}//在图表上显示滤波前后效果displayFilterResult(timeData, voltageData, filteredData, maxFrequency);return a.exec();
}QT滤波实现效果如下
