subplot(2,1,1); plot(fm, Zp_mag,'r*'); hold on; subplot(2,1,2); plot(fm, Zp_angle,'r*'); hold on; subplot(2,1,1); plot(fm, Zn_mag,'ro'); hold on; subplot(2,1,2); plot(fm, Zn_angle,'ro'); hold on;
时间: 2024-04-21 19:25:13 浏览: 99
这段代码使用了 `subplot` 函数创建了一个 2x1 的子图布局,并分别绘制了红色星号和红色圆圈的数据。
首先,在第一个子图中,使用 `subplot(2,1,1)` 将当前 axes 设置为第一个子图,然后使用 `plot` 函数绘制了 `fm` 和 `Zp_mag` 数据的散点图,颜色为红色星号。接着使用 `hold on` 保持绘图状态,以便在同一个子图上绘制其他数据。
然后,在第二个子图中,使用 `subplot(2,1,2)` 将当前 axes 设置为第二个子图,使用 `plot` 函数绘制了 `fm` 和 `Zp_angle` 数据的散点图,颜色同样为红色星号。同样地,使用 `hold on` 保持绘图状态。
接着,又回到第一个子图中,使用 `subplot(2,1,1)` 将当前 axes 设置为第一个子图,使用 `plot` 函数绘制了 `fm` 和 `Zn_mag` 数据的散点图,颜色为红色圆圈。再次使用 `hold on` 保持绘图状态。
最后,在第二个子图中,使用 `subplot(2,1,2)` 将当前 axes 设置为第二个子图,使用 `plot` 函数绘制了 `fm` 和 `Zn_angle` 数据的散点图,颜色同样为红色圆圈。最后再次使用 `hold on` 保持绘图状态。
这样,就在两个子图中分别绘制了红色星号和红色圆圈的数据。
相关问题
clear;clear;clc; fs = 8000; % 采样频率 t = 0:1/fs:1 ;% 采样时间 f = 2000; % 方波信号频率 x = square(2*pi*f*t); % 生成方波信号 N = length(x)-1; % 信号长度保证为正整数 % 计算基4FFT和频谱 X = fft(x, N/4); X_mag = abs(X); X_mag_norm = X_mag / max(X_mag); % 计算谐波失真度和谐波的归一化振幅 harmonics = 5; % 要计算的谐波次数 thd = 0; harmonic_amp_norm = zeros(1, harmonics); for k = 1:harmonics harmonic_amp_norm(k) = X_mag_norm(k*4+1); thd = thd + harmonic_amp_norm(k)^2; end thd = sqrt(thd) / harmonic_amp_norm(1) * 10; %打印失真度和归一化振幅图像 fprintf('前5个谐波分量的归一化振幅为:\n'); fprintf('%f\n', harmonic_amp_norm); fprintf('总谐波失真度:\n'); fprintf('%f\n',thd ); % 生成失真度和归一化振幅图像 figure(); subplot(211); stem(0:harmonics-1, harmonic_amp_norm); title('谐波振幅 (归一化)'); xlabel('谐波数'); ylabel('振幅'); subplot(212); plot(X_mag_norm); title(sprintf('总谐波失真度:THD = %.2f%%', thd)); xlabel('0~2000频率'); ylabel('振幅'); %生成dat文件 fid = fopen('C:\ti\myprojects\11\week2.dat','w'); fprintf(fid,'1651 2 0x8000 8000\n'); fprintf(fid,'%d\n',floor(x(1:N)*4095)); fclose(fid);通过生成的dat文件生成ccs中的C语言
您可以使用以下的 C 语言代码将生成的 `week2.dat` 文件读入到数组中,然后将数组中的数据写入到 DAC 寄存器中,以实现方波信号的输出:
```c
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <stdbool.h>
#include <math.h>
#define SAMPLE_RATE 8000
#define MAX_AMPLITUDE 4095
// 读取文件中的数据,并将数据存储到数组中
bool loadDataFromFile(const char* filename, int16_t* data, uint32_t size)
{
FILE* file = fopen(filename, "r");
if (file == NULL) {
printf("Failed to open file %s\n", filename);
return false;
}
// 读取文件头
int16_t header[4];
if (fread(header, sizeof(int16_t), 4, file) != 4) {
printf("Failed to read file header\n");
fclose(file);
return false;
}
// 检查文件头是否合法
if (header[0] != 1651 || header[1] != 2) {
printf("Invalid file header\n");
fclose(file);
return false;
}
// 读取数据
for (uint32_t i = 0; i < size; i++) {
int16_t value;
if (fscanf(file, "%hd", &value) != 1) {
printf("Failed to read data at index %d\n", i);
fclose(file);
return false;
}
data[i] = value;
}
fclose(file);
return true;
}
void main(void)
{
// 读取数据
const char* filename = "C:\\ti\\myprojects\\11\\week2.dat";
uint32_t dataSize = 8000;
int16_t data[dataSize];
if (!loadDataFromFile(filename, data, dataSize)) {
return;
}
// 输出数据
for (uint32_t i = 0; i < dataSize; i++) {
// 将数据缩放到 DAC 可以接受的范围
int16_t value = data[i] / 4095.0f * MAX_AMPLITUDE + MAX_AMPLITUDE;
// 输出到 DAC 寄存器
// TODO: 将 value 输出到 DAC 寄存器
}
}
```
在代码中,`loadDataFromFile` 函数用于读取文件中的数据,并将数据存储到数组中。`main` 函数中,首先读取 `week2.dat` 文件中的数据,然后对数据进行缩放,最后将数据输出到 DAC 寄存器中。由于不清楚您使用的具体的 DSP 芯片型号和开发板类型,因此无法提供具体的输出代码。您需要根据您的硬件和软件条件进行相应的修改,以实现数据的输出。
from scipy import signal import numpy as np import matplotlib.pyplot as plt plt.rcParams["font.family"] = 'Arial Unicode MS' original_sig = np.loadtxt("resources/unbalanced.txt") original_sig -= np.mean(original_sig) N = len(original_sig) pi = np.pi f2_jw = np.fft.fft(original_sig) f2_jw = np.fft.fftshift(f2_jw) jw_list = [complex(0, 1) * 2 * pi / N * item for item in np.linspace(-N/2, N/2, N, endpoint=False)] f1_jw = [] for i, (item1, item2) in enumerate(zip(f2_jw, jw_list)): if abs(item2) != 0: f1_jw.append(item1/item2) else: f1_jw.append(complex(0, 0)) f1_jw = np.array(f1_jw) * 1000 # m到mm的量纲转换 f1_jw = np.fft.ifftshift(f1_jw) vel_sig = np.fft.ifft(f1_jw).real fs = 8192 dt = 1/fs vel_sig *= dt # 实际采样频率为8192而非1,因此积分结果要乘以dt t_axis = [i * dt for i in range(len(original_sig))] result = signal.detrend(vel_sig) plt.figure(figsize=(12, 3)) plt.subplot(121) plt.plot(t_axis, vel_sig, label="频域积分计算得到的速度信号") plt.legend(loc="upper right") plt.subplot(122) plt.plot(t_axis, result, label="频域积分后去趋势得到的速度信号") plt.legend(loc="upper right") plt.show()将这段代码使用C语言进行编写,原始样本长度为512,为实数,在进行FFT处理之前,原始样本以一维数组表示复数,实部与虚部交替存在,虚部全部设置为0,FFT函数调用DSP库中的arm_cfft_f32函数
根据您提供的Python代码,我编写了相应的C语言代码,其中使用arm_cfft_f32函数进行FFT计算,并使用arm_rfft_f32函数进行实数FFT计算。同时,由于原始样本以一维数组表示复数,实部与虚部交替存在,虚部全部设置为0,因此在进行FFT计算时,需要将原始样本的虚部设置为0。具体实现如下:
```c
#include "arm_math.h"
#include <stdio.h>
#define N 512
int main()
{
// 读取原始样本
float original_sig[N];
FILE *fp = fopen("resources/unbalanced.txt", "r");
for (int i = 0; i < N; ++i) {
fscanf(fp, "%f", &original_sig[i]);
}
fclose(fp);
// 去趋势
arm_mean_f32(original_sig, N, &mean);
for (int i = 0; i < N; ++i) {
original_sig[i] -= mean;
}
// 使用一维数组表示复数,实部与虚部交替存在,虚部全部设置为0
float complex_sig[2 * N];
for (int i = 0; i < N; ++i) {
complex_sig[2 * i] = original_sig[i];
complex_sig[2 * i + 1] = 0;
}
// 进行FFT
arm_cfft_f32(&arm_cfft_sR_f32_len512, complex_sig, 0, 1);
arm_cmplx_mag_f32(complex_sig, original_sig, N);
// 进行频域积分
float jw_list[N];
float dt = 1.0 / 8192;
for (int i = 0; i < N; ++i) {
if (i < N / 2) {
jw_list[i] = 2 * ARM_PI / N * (i - N / 2);
} else {
jw_list[i] = 2 * ARM_PI / N * (i - N / 2);
}
}
for (int i = 0; i < N; ++i) {
if (i == 0) {
complex_sig[2 * i] = 0;
complex_sig[2 * i + 1] = 0;
} else {
complex_sig[2 * i] /= jw_list[i];
complex_sig[2 * i + 1] /= jw_list[i];
}
}
// 进行实数FFT计算
arm_rfft_instance_f32 rfft;
arm_rfft_init_f32(&rfft, N, 0, 1);
arm_rfft_f32(&rfft, complex_sig, original_sig);
// 转换为速度信号
for (int i = 0; i < N; ++i) {
original_sig[i] /= (2 * ARM_PI * (i == 0 ? 1 : jw_list[i]));
original_sig[i] *= 1000; // m到mm的量纲转换
original_sig[i] *= dt; // 实际采样频率为8192而非1,因此积分结果要乘以dt
}
// 去趋势
float mean;
arm_mean_f32(original_sig, N, &mean);
arm_offset_f32(original_sig, -mean, original_sig, N);
// 输出结果
for (int i = 0; i < N; ++i) {
printf("%f\n", original_sig[i]);
}
return 0;
}
```
需要注意的是,由于实数FFT计算得到的结果是对称的,因此在进行频域积分时,需要将虚部全部设置为0,并且jw_list需要进行特殊处理以保证积分结果的正确性。另外,为了方便计算,此处使用了arm_cmplx_mag_f32函数计算FFT的模值,而非手动计算模值。
阅读全文