解释 sig_fft = fft(signal,Nr)./Nr;
时间: 2024-05-24 09:10:47 浏览: 94
这行代码的作用是对长度为Nr的信号signal进行快速傅里叶变换(FFT),并将结果除以Nr得到归一化的FFT结果sig_fft。
其中,fft(signal,Nr)表示对信号signal进行FFT变换,Nr表示FFT的点数,即变换后得到的频域点数。./表示逐元素除法,将FFT结果中的每个元素都除以Nr,实现归一化。
相关问题
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
以下是使用C语言编写的代码:
```c
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define PI 3.14159265358979323846
void fft(double complex *x, int n) {
if (n == 1) return;
double complex xe[n/2], xo[n/2];
for (int i = 0; i < n/2; ++i) {
xe[i] = x[2*i];
xo[i] = x[2*i+1];
}
fft(xe, n/2);
fft(xo, n/2);
for (int i = 0; i < n/2; ++i) {
double complex t = cexp(-I * PI * i / (n/2)) * xo[i];
x[i] = xe[i] + t;
x[i+n/2] = xe[i] - t;
}
}
int main() {
double original_sig[512];
FILE *fp;
fp = fopen("resources/unbalanced.txt", "r");
for (int i = 0; i < 512; ++i) {
fscanf(fp, "%lf", &original_sig[i]);
original_sig[i] -= 0.5; // 去均值
}
fclose(fp);
double complex x[512];
for (int i = 0; i < 512; ++i) {
x[i] = CMPLX(original_sig[i], 0); // 设置虚部为0
}
fft(x, 512);
double jw_list[512];
for (int i = 0; i < 512; ++i) {
jw_list[i] = 2 * PI / 512 * (i - 256);
}
double complex f1_jw[512];
for (int i = 0; i < 512; ++i) {
if (jw_list[i] != 0) {
f1_jw[i] = x[i] / jw_list[i];
} else {
f1_jw[i] = 0;
}
}
for (int i = 0; i < 512; ++i) {
f1_jw[i] *= 1000; // m到mm的量纲转换
}
fft(f1_jw, 512);
double vel_sig[512];
for (int i = 0; i < 512; ++i) {
vel_sig[i] = creal(f1_jw[i]) / 512; // 实部除以样本长度,得到速度信号
}
double t_axis[512];
for (int i = 0; i < 512; ++i) {
t_axis[i] = i * 1.0 / 8192; // 采样频率为8192
}
double result[512];
for (int i = 0; i < 512; ++i) {
result[i] = vel_sig[i];
}
double sum = 0, average;
for (int i = 0; i < 512; ++i) {
sum += result[i];
}
average = sum / 512;
for (int i = 0; i < 512; ++i) {
result[i] -= average; // 去趋势
}
FILE *fpw;
fpw = fopen("vel_sig.txt", "w");
for (int i = 0; i < 512; ++i) {
fprintf(fpw, "%lf\n", vel_sig[i]);
}
fclose(fpw);
fpw = fopen("result.txt", "w");
for (int i = 0; i < 512; ++i) {
fprintf(fpw, "%lf\n", result[i]);
}
fclose(fpw);
return 0;
}
```
这段代码与Python代码实现的步骤基本相同,使用了C语言中的复数类型 `double complex` 和一些基本的数学函数。代码中还包含了读写文件的功能,将速度信号和去趋势后的速度信号分别写入了 `vel_sig.txt` 和 `result.txt` 文件中。
修改以下代码,plot向量长度不同% 读取输出文件数据 filename = 'D:\output_file.txt'; data = importdata(filename); Tp = 10e-6; B = 10e6; mu = B/Tp; % 采样率和时间轴 Fs = 400e6; Ts = 1/Fs; T_sample =Tp; N_sample = ceil(T_sample/Ts); t = (0:numel(data)-1) * Ts; Sig_ref = exp(1i*pi*mu*(t).^2); N_fft2 = length(data)+N_sample; PC_Sig_jam = fftshift( ifft(fft(data,N_fft2).*(conj(fft(Sig_ref,N_fft2)))) ); % 绘制脉冲压缩后的波形 figure; plot(t, PC_Sig_jam); xlabel('Time (s)'); ylabel('Amplitude'); title('Pulse Compression'); % 显示结果 disp('时域波形:'); disp(data); disp('脉冲压缩后的波形:'); disp(PC_Sig_jamd_data);
可以将原始数据和压缩后的数据都进行插值,使它们的长度相同,然后再进行绘图。
修改后的代码如下:
```matlab
% 读取输出文件数据
filename = 'D:\output_file.txt';
data = importdata(filename);
Tp = 10e-6;
B = 10e6;
mu = B/Tp;
% 采样率和时间轴
Fs = 400e6;
Ts = 1/Fs;
T_sample =Tp;
N_sample = ceil(T_sample/Ts);
t = (0:numel(data)-1) * Ts;
Sig_ref = exp(1i*pi*mu*(t).^2);
N_fft2 = length(data)+N_sample;
PC_Sig_jam = fftshift( ifft(fft(data,N_fft2).*(conj(fft(Sig_ref,N_fft2)))) );
% 插值
x = linspace(t(1), t(end), length(PC_Sig_jam));
y1 = interp1(t, data, x, 'linear');
y2 = interp1(t, PC_Sig_jam, x, 'linear');
% 绘制脉冲压缩后的波形
figure;
plot(x, y2);
xlabel('Time (s)');
ylabel('Amplitude');
title('Pulse Compression');
% 显示结果
disp('时域波形:');
disp(y1);
disp('脉冲压缩后的波形:');
disp(y2);
```
这样就可以避免向量长度不同的问题。
阅读全文