心电信号去除50hz工频干扰

% 读取文本文件中的心电信号数据
filename = 'ECG_X1.txt';
data = load(filename);

% 绘制心电信号的时域波形
figure;
plot(data);
axis("tight");
title('心电信号时域波形');
xlabel('采样点');
ylabel('幅值');

% 计算心电信号的FFT并绘制频谱图
N = length(data); % 信号长度
Fs = 150; % 采样率（假设为1000Hz，根据实际情况修改）
frequencies = Fs*(0:(N/2))/N; % 计算频率轴
fft_data = fft(data); % 计算FFT
fft_data_abs = abs(fft_data/N); % 计算幅值
fft_data_abs_half = fft_data_abs(1:N/2+1); % 取FFT结果的前一半

figure;
plot(frequencies, fft_data_abs_half);
title('心电信号频谱图');
xlabel('频率 (Hz)');
ylabel('幅值');

% 找出频谱图中的主要频率成分
[max_amp, max_index] = max(fft_data_abs_half);
dominant_freq = frequencies(max_index);
fprintf('主要频率成分为 %.2f Hz\n', dominant_freq);



%%
t = (0:length(data)-1) / Fs; % 生成时间向量
noise_amplitude = 20; % 干扰信号的振幅
noise = noise_amplitude * sin(2*pi*50*t); % 生成50Hz的工频干扰信号
data_with_noise = data' + noise; % 添加干扰信号

% 绘制添加干扰信号后的心电信号时域波形
figure;
plot(t, data_with_noise);
title('添加50Hz工频干扰后的心电信号时域波形');
xlabel('时间 (s)');
ylabel('幅值');

% 计算心电信号加干扰的FFT并绘制频谱图
N = length(data_with_noise); % 信号长度
frequencies = Fs*(0:(N/2))/N; % 计算频率轴
fft_data_with_noise = fft(data_with_noise); % 计算FFT
fft_data_with_noise_abs = abs(fft_data_with_noise/N); % 计算幅值
fft_data_with_noise_abs_half = fft_data_with_noise_abs(1:N/2+1); % 取FFT结果的前一半

% 绘制心电信号加干扰的频谱图
figure;
plot(frequencies, fft_data_with_noise_abs_half);
title('心电信号加50Hz工频干扰的频谱图');
xlabel('频率 (Hz)');
ylabel('幅值');

% 找出频谱图中的主要频率成分
[max_amp, max_index] = max(fft_data_with_noise_abs_half);
dominant_freq = frequencies(max_index);
fprintf('主要频率成分为 %.2f Hz\n', dominant_freq);

%% FIR 
Fs2=Fs/2;
fp = 35; fs = 45;            % 通带和阻带频率
wp = fp*pi/Fs2; ws = fs*pi/Fs2;  % 通带和阻带归一化角频率
deltaw= ws - wp;                 % 过渡带宽Δω的计算

N = ceil(6.6*pi/ deltaw);        % 计算所需的滤波器阶数N
N = N + mod(N,2);                % 保证滤波器系数为奇数

wind = (hamming(N+1))';          % 海明窗计算
Wn=(fp+fs)/2/Fs2;                % 计算截止频率

b = fir1(N,Wn,wind);             % 用fir1函数设计滤波器
figure();freqz(b,1);             % 滤波器画图

data_fir_out= filter(b,1,data_with_noise);
figure;
plot(t, out);
title('FIR滤波去除50Hz工频干扰后的心电信号时域波形');
xlabel('时间 (s)');
ylabel('幅值');

% 计算心电信号去除干扰的FFT并绘制频谱图
N = length(data_fir_out); % 信号长度
frequencies = Fs*(0:(N/2))/N; % 计算频率轴
fft_data_fir_out = fft(data_fir_out); % 计算FFT
fft_data_fir_out_abs = abs(fft_data_fir_out/N); % 计算幅值
fft_data_fir_out_abs_half = fft_data_fir_out_abs(1:N/2+1); % 取FFT结果的前一半

% 绘制FIR滤波后的频谱图
figure;
plot(frequencies, fft_data_fir_out_abs_half);
title('心电信号通过FIR滤波后频谱图');
xlabel('频率 (Hz)');
ylabel('幅值');

%% LMS自适应滤波