【matlab】读取rhd文件格式，intan脑电采集开源代码
MATLAB 脚本文档：read_Intan_RHD2000_file
概述
此 MATLAB 脚本 read_Intan_RHD2000_file 旨在读取由 Intan Technologies RHD2000 系统生成的数据文件。它解析数据并将其置于 MATLAB 工作空间中，以便进一步分析。脚本兼容 Intan USB 接口板和 Intan 录音控制器。
先决条件
- 系统上安装了 MATLAB。
- RHD2000 数据文件（
*.rhd）在文件系统中可访问。
使用方法
-
清除工作空间：在运行脚本之前，建议清除 MATLAB 工作空间，以避免与现有变量冲突。
clear -
运行脚本：在 MATLAB 中执行脚本。
read_Intan_RHD2000_file -
选择数据文件：文件选择对话框将提示您选择一个 RHD2000 数据文件（
*.rhd）。选择文件并点击“打开”。 -
查看变量：脚本完成后，您可以使用
whos命令查看 MATLAB 工作空间中的可用变量。whos -
访问数据：您可以直接从工作空间访问解析后的数据。例如，要访问第一个放大器通道数据：
amplifier_data(1, :)
主要特性
- 自动文件检测：脚本自动检测并读取文件，基于“魔术数字”确保它是 RHD2000 数据文件。
- 版本兼容性：支持不同版本的 RHD2000 数据格式。
- 数据解析：提取各种数据类型，包括放大器通道、辅助输入通道、供电电压通道以及板载 ADC/DAC 通道。
- 滤波器应用：如果录音时启用了此类滤波器，则可以选择对放大器数据应用陷波滤波器。
- 数据缩放：将数据缩放到适当的单位（例如，放大器数据为微伏，辅助输入为伏特）。
数据结构
- 放大器通道：
amplifier_channels,amplifier_data - 辅助输入通道：
aux_input_channels,aux_input_data - 供电电压通道：
supply_voltage_channels,supply_voltage_data - 板载 ADC 通道：
board_adc_channels,board_adc_data - 板载数字输入/输出通道：
board_dig_in_channels,board_dig_out_channels,board_dig_in_data,board_dig_out_data - 时间戳：
t_amplifier,t_aux_input,t_supply_voltage,t_board_adc,t_dig,t_temp_sensor
附加实用工具
- 陷波滤波器函数：
notch_filter- 对数据应用陷波滤波器。 - 复数化函数：
plural- 根据计数为单词添加 ‘s’ 的实用工具函数。 - QString 读取器：
fread_QString- 从文件中读取 Qt 风格的 QString。
输出
- 脚本将进度和结果输出到 MATLAB 命令窗口，包括数据读取、解析的状态以及遇到的任何警告或错误。
故障排除
- 确保 RHD2000 数据文件未损坏且与脚本兼容。
- 如果遇到错误，请检查 MATLAB 版本兼容性。
- 验证文件路径是否正确且文件可访问。
结论
read_Intan_RHD2000_file 脚本是导入和分析 RHD2000 数据的全面工具。它为处理数据提供了结构化的方法，使研究人员和工程师更容易进一步分析和可视化数据。
代码
function read_Intan_RHD2000_file
% read_Intan_RHD2000_file
%
% Version 3.0, 8 February 2021
%
% Reads Intan Technologies RHD data file generated by Intan USB interface
% board or Intan Recording Controller. Data are parsed and placed into
% variables that appear in the base MATLAB workspace. Therefore, it is
% recommended to execute a 'clear' command before running this program to
% clear all other variables from the base workspace.
%
% Example:
% >> clear
% >> read_Intan_RHD2000_file
% >> whos
% >> amplifier_channels(1)
% >> plot(t_amplifier, amplifier_data(1,:))
[file, path, filterindex] = ...
uigetfile('*.rhd', 'Select an RHD2000 Data File', 'MultiSelect', 'off');
if (file == 0)
return;
end
tic;
filename = [path,file];
fid = fopen(filename, 'r');
s = dir(filename);
filesize = s.bytes;
% Check 'magic number' at beginning of file to make sure this is an Intan
% Technologies RHD2000 data file.
magic_number = fread(fid, 1, 'uint32');
if magic_number ~= hex2dec('c6912702')
error('Unrecognized file type.');
end
% Read version number.
data_file_main_version_number = fread(fid, 1, 'int16');
data_file_secondary_version_number = fread(fid, 1, 'int16');
fprintf(1, '\n');
fprintf(1, 'Reading Intan Technologies RHD2000 Data File, Version %d.%d\n', ...
data_file_main_version_number, data_file_secondary_version_number);
fprintf(1, '\n');
if (data_file_main_version_number == 1)
num_samples_per_data_block = 60;
else
num_samples_per_data_block = 128;
end
% Read information of sampling rate and amplifier frequency settings.
sample_rate = fread(fid, 1, 'single');
dsp_enabled = fread(fid, 1, 'int16');
actual_dsp_cutoff_frequency = fread(fid, 1, 'single');
actual_lower_bandwidth = fread(fid, 1, 'single');
actual_upper_bandwidth = fread(fid, 1, 'single');
desired_dsp_cutoff_frequency = fread(fid, 1, 'single');
desired_lower_bandwidth = fread(fid, 1, 'single');
desired_upper_bandwidth = fread(fid, 1, 'single');
% This tells us if a software 50/60 Hz notch filter was enabled during
% the data acquisition.
notch_filter_mode = fread(fid, 1, 'int16');
notch_filter_frequency = 0;
if (notch_filter_mode == 1)
notch_filter_frequency = 50;
elseif (notch_filter_mode == 2)
notch_filter_frequency = 60;
end
desired_impedance_test_frequency = fread(fid, 1, 'single');
actual_impedance_test_frequency = fread(fid, 1, 'single');
% Place notes in data strucure
notes = struct( ...
'note1', fread_QString(fid), ...
'note2', fread_QString(fid), ...
'note3', fread_QString(fid) );
% If data file is from GUI v1.1 or later, see if temperature sensor data
% was saved.
num_temp_sensor_channels = 0;
if ((data_file_main_version_number == 1 && data_file_secondary_version_number >= 1) ...
|| (data_file_main_version_number > 1))
num_temp_sensor_channels = fread(fid, 1, 'int16');
end
% If data file is from GUI v1.3 or later, load board mode.
board_mode = 0;
if ((data_file_main_version_number == 1 && data_file_secondary_version_number >= 3) ...
|| (data_file_main_version_number > 1))
board_mode = fread(fid, 1, 'int16');
end
% If data file is from v2.0 or later (Intan Recording Controller),
% load name of digital reference channel.
if (data_file_main_version_number > 1)
reference_channel = fread_QString(fid);
end
% Place frequency-related information in data structure.
frequency_parameters = struct( ...
'amplifier_sample_rate', sample_rate, ...
'aux_input_sample_rate', sample_rate / 4, ...
'supply_voltage_sample_rate', sample_rate / num_samples_per_data_block, ...
'board_adc_sample_rate', sample_rate, ...
'board_dig_in_sample_rate', sample_rate, ...
'desired_dsp_cutoff_frequency', desired_dsp_cutoff_frequency, ...
'actual_dsp_cutoff_frequency', actual_dsp_cutoff_frequency, ...
'dsp_enabled', dsp_enabled, ...
'desired_lower_bandwidth', desired_lower_bandwidth, ...
'actual_lower_bandwidth', actual_lower_bandwidth, ...
'desired_upper_bandwidth', desired_upper_bandwidth, ...
'actual_upper_bandwidth', actual_upper_bandwidth, ...
'notch_filter_frequency', notch_filter_frequency, ...
'desired_impedance_test_frequency', desired_impedance_test_frequency, ...
'actual_impedance_test_frequency', actual_impedance_test_frequency );
% Define data structure for spike trigger settings.
spike_trigger_struct = struct( ...
'voltage_trigger_mode', {}, ...
'voltage_threshold', {}, ...
'digital_trigger_channel', {}, ...
'digital_edge_polarity', {} );
new_trigger_channel = struct(spike_trigger_struct);
spike_triggers = struct(spike_trigger_struct);
% Define data structure for data channels.
channel_struct = struct( ...
'native_channel_name', {}, ...
'custom_channel_name', {}, ...
'native_order', {}, ...
'custom_order', {}, ...
'board_stream', {}, ...
'chip_channel', {}, ...
'port_name', {}, ...
'port_prefix', {}, ...
'port_number', {}, ...
'electrode_impedance_magnitude', {}, ...
'electrode_impedance_phase', {} );
new_channel = struct(channel_struct);
% Create structure arrays for each type of data channel.
amplifier_channels = struct(channel_struct);
aux_input_channels = struct(channel_struct);
supply_voltage_channels = struct(channel_struct);
board_adc_channels = struct(channel_struct);
board_dig_in_channels = struct(channel_struct);
board_dig_out_channels = struct(channel_struct);
amplifier_index = 1;
aux_input_index = 1;
supply_voltage_index = 1;
board_adc_index = 1;
board_dig_in_index = 1;
board_dig_out_index = 1;
% Read signal summary from data file header.
number_of_signal_groups = fread(fid, 1, 'int16');
for signal_group = 1:number_of_signal_groups
signal_group_name = fread_QString(fid);
signal_group_prefix = fread_QString(fid);
signal_group_enabled = fread(fid, 1, 'int16');
signal_group_num_channels = fread(fid, 1, 'int16');
signal_group_num_amp_channels = fread(fid, 1, 'int16');
if (signal_group_num_channels > 0 && signal_group_enabled > 0)
new_channel(1).port_name = signal_group_name;
new_channel(1).port_prefix = signal_group_prefix;
new_channel(1).port_number = signal_group;
for signal_channel = 1:signal_group_num_channels
new_channel(1).native_channel_name = fread_QString(fid);
new_channel(1).custom_channel_name = fread_QString(fid);
new_channel(1).native_order = fread(fid, 1, 'int16');
new_channel(1).custom_order = fread(fid, 1, 'int16');
signal_type = fread(fid, 1, 'int16');
channel_enabled = fread(fid, 1, 'int16');
new_channel(1).chip_channel = fread(fid, 1, 'int16');
new_channel(1).board_stream = fread(fid, 1, 'int16');
new_trigger_channel(1).voltage_trigger_mode = fread(fid, 1, 'int16');
new_trigger_channel(1).voltage_threshold = fread(fid, 1, 'int16');
new_trigger_channel(1).digital_trigger_channel = fread(fid, 1, 'int16');
new_trigger_channel(1).digital_edge_polarity = fread(fid, 1, 'int16');
new_channel(1).electrode_impedance_magnitude = fread(fid, 1, 'single');
new_channel(1).electrode_impedance_phase = fread(fid, 1, 'single');
if (channel_enabled)
switch (signal_type)
case 0
amplifier_channels(amplifier_index) = new_channel;
spike_triggers(amplifier_index) = new_trigger_channel;
amplifier_index = amplifier_index + 1;
case 1
aux_input_channels(aux_input_index) = new_channel;
aux_input_index = aux_input_index + 1;
case 2
supply_voltage_channels(supply_voltage_index) = new_channel;
supply_voltage_index = supply_voltage_index + 1;
case 3
board_adc_channels(board_adc_index) = new_channel;
board_adc_index = board_adc_index + 1;
case 4
board_dig_in_channels(board_dig_in_index) = new_channel;
board_dig_in_index = board_dig_in_index + 1;
case 5
board_dig_out_channels(board_dig_out_index) = new_channel;
board_dig_out_index = board_dig_out_index + 1;
otherwise
error('Unknown channel type');
end
end
end
end
end
% Summarize contents of data file.
num_amplifier_channels = amplifier_index - 1;
num_aux_input_channels = aux_input_index - 1;
num_supply_voltage_channels = supply_voltage_index - 1;
num_board_adc_channels = board_adc_index - 1;
num_board_dig_in_channels = board_dig_in_index - 1;
num_board_dig_out_channels = board_dig_out_index - 1;
fprintf(1, 'Found %d amplifier channel%s.\n', ...
num_amplifier_channels, plural(num_amplifier_channels));
fprintf(1, 'Found %d auxiliary input channel%s.\n', ...
num_aux_input_channels, plural(num_aux_input_channels));
fprintf(1, 'Found %d supply voltage channel%s.\n', ...
num_supply_voltage_channels, plural(num_supply_voltage_channels));
fprintf(1, 'Found %d board ADC channel%s.\n', ...
num_board_adc_channels, plural(num_board_adc_channels));
fprintf(1, 'Found %d board digital input channel%s.\n', ...
num_board_dig_in_channels, plural(num_board_dig_in_channels));
fprintf(1, 'Found %d board digital output channel%s.\n', ...
num_board_dig_out_channels, plural(num_board_dig_out_channels));
fprintf(1, 'Found %d temperature sensor channel%s.\n', ...
num_temp_sensor_channels, plural(num_temp_sensor_channels));
fprintf(1, '\n');
% Determine how many samples the data file contains.
% Each data block contains num_samples_per_data_block amplifier samples.
bytes_per_block = num_samples_per_data_block * 4; % timestamp data
bytes_per_block = bytes_per_block + num_samples_per_data_block * 2 * num_amplifier_channels;
% Auxiliary inputs are sampled 4x slower than amplifiers
bytes_per_block = bytes_per_block + (num_samples_per_data_block / 4) * 2 * num_aux_input_channels;
% Supply voltage is sampled once per data block
bytes_per_block = bytes_per_block + 1 * 2 * num_supply_voltage_channels;
% Board analog inputs are sampled at same rate as amplifiers
bytes_per_block = bytes_per_block + num_samples_per_data_block * 2 * num_board_adc_channels;
% Board digital inputs are sampled at same rate as amplifiers
if (num_board_dig_in_channels > 0)
bytes_per_block = bytes_per_block + num_samples_per_data_block * 2;
end
% Board digital outputs are sampled at same rate as amplifiers
if (num_board_dig_out_channels > 0)
bytes_per_block = bytes_per_block + num_samples_per_data_block * 2;
end
% Temp sensor is sampled once per data block
if (num_temp_sensor_channels > 0)
bytes_per_block = bytes_per_block + 1 * 2 * num_temp_sensor_channels;
end
% How many data blocks remain in this file?
data_present = 0;
bytes_remaining = filesize - ftell(fid);
if (bytes_remaining > 0)
data_present = 1;
end
num_data_blocks = bytes_remaining / bytes_per_block;
num_amplifier_samples = num_samples_per_data_block * num_data_blocks;
num_aux_input_samples = (num_samples_per_data_block / 4) * num_data_blocks;
num_supply_voltage_samples = 1 * num_data_blocks;
num_board_adc_samples = num_samples_per_data_block * num_data_blocks;
num_board_dig_in_samples = num_samples_per_data_block * num_data_blocks;
num_board_dig_out_samples = num_samples_per_data_block * num_data_blocks;
record_time = num_amplifier_samples / sample_rate;
if (data_present)
fprintf(1, 'File contains %0.3f seconds of data. Amplifiers were sampled at %0.2f kS/s.\n', ...
record_time, sample_rate / 1000);
fprintf(1, '\n');
else
fprintf(1, 'Header file contains no data. Amplifiers were sampled at %0.2f kS/s.\n', ...
sample_rate / 1000);
fprintf(1, '\n');
end
if (data_present)
% Pre-allocate memory for data.
fprintf(1, 'Allocating memory for data...\n');
t_amplifier = zeros(1, num_amplifier_samples);
amplifier_data = zeros(num_amplifier_channels, num_amplifier_samples);
aux_input_data = zeros(num_aux_input_channels, num_aux_input_samples);
supply_voltage_data = zeros(num_supply_voltage_channels, num_supply_voltage_samples);
temp_sensor_data = zeros(num_temp_sensor_channels, num_supply_voltage_samples);
board_adc_data = zeros(num_board_adc_channels, num_board_adc_samples);
board_dig_in_data = zeros(num_board_dig_in_channels, num_board_dig_in_samples);
board_dig_in_raw = zeros(1, num_board_dig_in_samples);
board_dig_out_data = zeros(num_board_dig_out_channels, num_board_dig_out_samples);
board_dig_out_raw = zeros(1, num_board_dig_out_samples);
% Read sampled data from file.
fprintf(1, 'Reading data from file...\n');
amplifier_index = 1;
aux_input_index = 1;
supply_voltage_index = 1;
board_adc_index = 1;
board_dig_in_index = 1;
board_dig_out_index = 1;
print_increment = 10;
percent_done = print_increment;
for i=1:num_data_blocks
% In version 1.2, we moved from saving timestamps as unsigned
% integeters to signed integers to accomidate negative (adjusted)
% timestamps for pretrigger data.
if ((data_file_main_version_number == 1 && data_file_secondary_version_number >= 2) ...
|| (data_file_main_version_number > 1))
t_amplifier(amplifier_index:(amplifier_index + num_samples_per_data_block - 1)) = fread(fid, num_samples_per_data_block, 'int32');
else
t_amplifier(amplifier_index:(amplifier_index + num_samples_per_data_block - 1)) = fread(fid, num_samples_per_data_block, 'uint32');
end
if (num_amplifier_channels > 0)
amplifier_data(:, amplifier_index:(amplifier_index + num_samples_per_data_block - 1)) = fread(fid, [num_samples_per_data_block, num_amplifier_channels], 'uint16')';
end
if (num_aux_input_channels > 0)
aux_input_data(:, aux_input_index:(aux_input_index + (num_samples_per_data_block / 4) - 1)) = fread(fid, [(num_samples_per_data_block / 4), num_aux_input_channels], 'uint16')';
end
if (num_supply_voltage_channels > 0)
supply_voltage_data(:, supply_voltage_index) = fread(fid, [1, num_supply_voltage_channels], 'uint16')';
end
if (num_temp_sensor_channels > 0)
temp_sensor_data(:, supply_voltage_index) = fread(fid, [1, num_temp_sensor_channels], 'int16')';
end
if (num_board_adc_channels > 0)
board_adc_data(:, board_adc_index:(board_adc_index + num_samples_per_data_block - 1)) = fread(fid, [num_samples_per_data_block, num_board_adc_channels], 'uint16')';
end
if (num_board_dig_in_channels > 0)
board_dig_in_raw(board_dig_in_index:(board_dig_in_index + num_samples_per_data_block - 1)) = fread(fid, num_samples_per_data_block, 'uint16');
end
if (num_board_dig_out_channels > 0)
board_dig_out_raw(board_dig_out_index:(board_dig_out_index + num_samples_per_data_block - 1)) = fread(fid, num_samples_per_data_block, 'uint16');
end
amplifier_index = amplifier_index + num_samples_per_data_block;
aux_input_index = aux_input_index + (num_samples_per_data_block / 4);
supply_voltage_index = supply_voltage_index + 1;
board_adc_index = board_adc_index + num_samples_per_data_block;
board_dig_in_index = board_dig_in_index + num_samples_per_data_block;
board_dig_out_index = board_dig_out_index + num_samples_per_data_block;
fraction_done = 100 * (i / num_data_blocks);
if (fraction_done >= percent_done)
fprintf(1, '%d%% done...\n', percent_done);
percent_done = percent_done + print_increment;
end
end
% Make sure we have read exactly the right amount of data.
bytes_remaining = filesize - ftell(fid);
if (bytes_remaining ~= 0)
%error('Error: End of file not reached.');
end
end
% Close data file.
fclose(fid);
if (data_present)
fprintf(1, 'Parsing data...\n');
% Extract digital input channels to separate variables.
for i=1:num_board_dig_in_channels
mask = 2^(board_dig_in_channels(i).native_order) * ones(size(board_dig_in_raw));
board_dig_in_data(i, :) = (bitand(board_dig_in_raw, mask) > 0);
end
for i=1:num_board_dig_out_channels
mask = 2^(board_dig_out_channels(i).native_order) * ones(size(board_dig_out_raw));
board_dig_out_data(i, :) = (bitand(board_dig_out_raw, mask) > 0);
end
% Scale voltage levels appropriately.
amplifier_data = 0.195 * (amplifier_data - 32768); % units = microvolts
aux_input_data = 37.4e-6 * aux_input_data; % units = volts
supply_voltage_data = 74.8e-6 * supply_voltage_data; % units = volts
if (board_mode == 1)
board_adc_data = 152.59e-6 * (board_adc_data - 32768); % units = volts
elseif (board_mode == 13) % Intan Recording Controller
board_adc_data = 312.5e-6 * (board_adc_data - 32768); % units = volts
else
board_adc_data = 50.354e-6 * board_adc_data; % units = volts
end
temp_sensor_data = temp_sensor_data / 100; % units = deg C
% Check for gaps in timestamps.
num_gaps = sum(diff(t_amplifier) ~= 1);
if (num_gaps == 0)
fprintf(1, 'No missing timestamps in data.\n');
else
fprintf(1, 'Warning: %d gaps in timestamp data found. Time scale will not be uniform!\n', ...
num_gaps);
end
% Scale time steps (units = seconds).
t_amplifier = t_amplifier / sample_rate;
t_aux_input = t_amplifier(1:4:end);
t_supply_voltage = t_amplifier(1:num_samples_per_data_block:end);
t_board_adc = t_amplifier;
t_dig = t_amplifier;
t_temp_sensor = t_supply_voltage;
% If the software notch filter was selected during the recording, apply the
% same notch filter to amplifier data here. But don't do this for v3.0+
% files (from Intan RHX software) because RHX saves notch-filtered data.
if (notch_filter_frequency > 0 && data_file_main_version_number < 3)
fprintf(1, 'Applying notch filter...\n');
print_increment = 10;
percent_done = print_increment;
for i=1:num_amplifier_channels
amplifier_data(i,:) = ...
notch_filter(amplifier_data(i,:), sample_rate, notch_filter_frequency, 10);
fraction_done = 100 * (i / num_amplifier_channels);
if (fraction_done >= percent_done)
fprintf(1, '%d%% done...\n', percent_done);
percent_done = percent_done + print_increment;
end
end
end
end
% Move variables to base workspace.
% new for version 2.01: move filename info to base workspace
filename = file;
move_to_base_workspace(filename);
move_to_base_workspace(path);
move_to_base_workspace(notes);
move_to_base_workspace(frequency_parameters);
if (data_file_main_version_number > 1)
move_to_base_workspace(reference_channel);
end
if (num_amplifier_channels > 0)
move_to_base_workspace(amplifier_channels);
if (data_present)
move_to_base_workspace(amplifier_data);
move_to_base_workspace(t_amplifier);
end
move_to_base_workspace(spike_triggers);
end
if (num_aux_input_channels > 0)
move_to_base_workspace(aux_input_channels);
if (data_present)
move_to_base_workspace(aux_input_data);
move_to_base_workspace(t_aux_input);
end
end
if (num_supply_voltage_channels > 0)
move_to_base_workspace(supply_voltage_channels);
if (data_present)
move_to_base_workspace(supply_voltage_data);
move_to_base_workspace(t_supply_voltage);
end
end
if (num_board_adc_channels > 0)
move_to_base_workspace(board_adc_channels);
if (data_present)
move_to_base_workspace(board_adc_data);
move_to_base_workspace(t_board_adc);
end
end
if (num_board_dig_in_channels > 0)
move_to_base_workspace(board_dig_in_channels);
if (data_present)
move_to_base_workspace(board_dig_in_data);
move_to_base_workspace(t_dig);
end
end
if (num_board_dig_out_channels > 0)
move_to_base_workspace(board_dig_out_channels);
if (data_present)
move_to_base_workspace(board_dig_out_data);
move_to_base_workspace(t_dig);
end
end
if (num_temp_sensor_channels > 0)
if (data_present)
move_to_base_workspace(temp_sensor_data);
move_to_base_workspace(t_temp_sensor);
end
end
fprintf(1, 'Done! Elapsed time: %0.1f seconds\n', toc);
if (data_present)
fprintf(1, 'Extracted data are now available in the MATLAB workspace.\n');
else
fprintf(1, 'Extracted waveform information is now available in the MATLAB workspace.\n');
end
fprintf(1, 'Type ''whos'' to see variables.\n');
fprintf(1, '\n');
return
function a = fread_QString(fid)
% a = read_QString(fid)
%
% Read Qt style QString. The first 32-bit unsigned number indicates
% the length of the string (in bytes). If this number equals 0xFFFFFFFF,
% the string is null.
a = '';
length = fread(fid, 1, 'uint32');
if length == hex2num('ffffffff')
return;
end
% convert length from bytes to 16-bit Unicode words
length = length / 2;
for i=1:length
a(i) = fread(fid, 1, 'uint16');
end
return
function s = plural(n)
% s = plural(n)
%
% Utility function to optionally plurailze words based on the value
% of n.
if (n == 1)
s = '';
else
s = 's';
end
return
function out = notch_filter(in, fSample, fNotch, Bandwidth)
% out = notch_filter(in, fSample, fNotch, Bandwidth)
%
% Implements a notch filter (e.g., for 50 or 60 Hz) on vector 'in'.
% fSample = sample rate of data (in Hz or Samples/sec)
% fNotch = filter notch frequency (in Hz)
% Bandwidth = notch 3-dB bandwidth (in Hz). A bandwidth of 10 Hz is
% recommended for 50 or 60 Hz notch filters; narrower bandwidths lead to
% poor time-domain properties with an extended ringing response to
% transient disturbances.
%
% Example: If neural data was sampled at 30 kSamples/sec
% and you wish to implement a 60 Hz notch filter:
%
% out = notch_filter(in, 30000, 60, 10);
tstep = 1/fSample;
Fc = fNotch*tstep;
L = length(in);
% Calculate IIR filter parameters
d = exp(-2*pi*(Bandwidth/2)*tstep);
b = (1 + d*d)*cos(2*pi*Fc);
a0 = 1;
a1 = -b;
a2 = d*d;
a = (1 + d*d)/2;
b0 = 1;
b1 = -2*cos(2*pi*Fc);
b2 = 1;
out = zeros(size(in));
out(1) = in(1);
out(2) = in(2);
% (If filtering a continuous data stream, change out(1) and out(2) to the
% previous final two values of out.)
% Run filter
for i=3:L
out(i) = (a*b2*in(i-2) + a*b1*in(i-1) + a*b0*in(i) - a2*out(i-2) - a1*out(i-1))/a0;
end
return
function move_to_base_workspace(variable)
% move_to_base_workspace(variable)
%
% Move variable from function workspace to base MATLAB workspace so
% user will have access to it after the program ends.
variable_name = inputname(1);
assignin('base', variable_name, variable);
return;