Digital FIR Filter Design

cancelling-effects.sci

close() ; close(); clear()

function [H] = compute_irc(signal)
    // spectrum = real(fft(signal));
    H = signal .* window('kr', length(signal), 16);
endfunction

function [out] = compute_inv_irc(signal)
    filter_fft = fft(signal)

    // N = length(filter_fft)
    // offsets = exp(-(%i * 2 * %pi * (-0.5 + (0:N-1) / N) * floor(N / 2)))
    
    H = 1 ./ filter_fft // .* offsets ; // fftshift(filter_fft) // .* offsets;
    H(1, find(isinf(H))) = 0.;
    // H = conj(filter_fft)./abs(filter_fft)^2

    out = ifft(H) // .* window('kr', length(H) , 8);
    out_len = length(out)
    mid = (out_len - modulo(out_len, 2)) / 2;
    out = cat(2, out(mid:out_len), out(1:mid-1));

    out = out .* window('kr', length(out), 16);
endfunction

function plot_spectrum(figure, signal, Fs)
    signal_len = length(signal) 
    signal_freq = (1:signal_len) / signal_len * Fs;
    scf(figure)

    magnitudes = fft(signal)

    plot2d(signal_freq, abs(magnitudes))
    title('frequencies')
endfunction

function plot_phase_response(figure, spectrum)
    [db,phi] = dbphi(spectrum)

    scf(figure)
    plot2d(phi)
endfunction

function [H] = apply_filter(signal, filter, do_fft)
    if ~exists("do_fft", "local") then
      do_fft = 1
    end

    if do_fft
        disp('computing conv using our fft')
        H = fft_convolution(signal, filter)
    else
        disp('computing conv using built-in')
        H = conv(signal, filter);
    end
endfunction

function [conv] = fft_convolution(sig, filt)
    tmp_size = 2^ceil(log2(size(sig, 'c')))

    filt_f = resize_matrix(filt, -1, tmp_size)
    sig_f = resize_matrix(sig, -1, tmp_size)

    out_size = size(sig, 'c') + size(filt, 'c')

    conv = resize_matrix(ifft(fft(filt_f).*fft(sig_f)), -1, out_size)
endfunction


[signal, Fs, _] = wavread('1a_marble_hall.wav'); // 'irc2.wav');
signal = signal(1, :);

// plot2d(signal)

subplot(331)
plot_spectrum(0, signal, Fs)

// irc = compute_irc(signal);
irc = signal;


subplot(332)
plot2d(irc / norm(irc, %inf))
title("irc")

[input, Fs, _] = wavread('7cef8230.wav');
input = input(1, :); // 1:length(input) / 64);

subplot(333)
plot2d(input)
title('orig signal')

// subplot(224)
convolved = apply_filter(input, irc, do_fft=1);
// plot2d(convolved)
// title('convolved')

subplot(334)
plot2d(convolved/norm(convolved, %inf))
title('convolved, normalized')

wavwrite(convolved, 'convolved.wav')
savewave('convolved.wav', convolved/norm(convolved, %inf), Fs)

inverse = compute_inv_irc(irc);

subplot(337)
plot_phase_response(0, fft(irc));
title('phase irc')
subplot(338)
plot_phase_response(0, fft(inverse));
title('phase inverse irc')

// close()

scf(1)
subplot(331)
plot2d(convolved)
title('convolved, normalized')

subplot(332)
original = apply_filter(convolved, inverse, do_fft=1);
plot2d(original/norm(original, %inf))
title('convolved with inverse irc')

savewave('restored.wav', original/norm(original, %inf), Fs)

subplot(333)
plot_spectrum(1, inverse, Fs)
//  plot2d(real((fft(inverse))))

subplot(334)
plot2d(inverse/norm(inverse, %inf))
title('inverse irc')

subplot(335)
plot2d(apply_filter(inverse, irc, do_fft=1))
title('approx of delta function')

fir-filter-design.sci

clear();

function H = ideal_lowpass(N, cutoff, stop_value)
    N = (N - modulo(N, 2)) / 2;
    cutoff = floor(2 * N * cutoff);
    H = ones(1, N) * stop_value
    H(1,1:cutoff) = 1.;
    // need to make N odd
    H = [1. H flipdim(H, 2)] // <---- line 14
endfunction

function H = ideal_highpass(N, low, high, stop_value)
    N = (N - modulo(N, 2)) / 2;
    cutoff_low = floor(2 * N * low);
    cutoff_high = floor(2 * N * high);

    H = ones(1, N) * stop_value;
    // H(1,1:cutoff_low) = 1.;

    write(%io(2), cutoff_high)
    write(%io(2), N)

    if (2 > cutoff_high)
        H(1,2:N) = 1
    else
        H(1,cutoff_high:N) = 1.;
    end

    H = [1. H flipdim(H, 2)] // <---- line 14
endfunction


function draw_freq_w_noise(signal, frequencies)
    magnitude = abs(fft(signal));
    plot2d(frequencies, magnitude, logflag='nl')
    xlabel("Freq, n", 'fontsize', 2)
    ylabel("Amplitude", 'fontsize', 2)
    title("Signal with noise in freq domain", 'fontsize', 3)
endfunction

function draw_ir_freq(signal, frequencies)
    plot2d(frequencies, signal, logflag='nn')
    xlabel("Freq, n", 'fontsize', 2)
    ylabel("Amplitude", 'fontsize', 2)
    title("IR in freq domain", 'fontsize', 3)
endfunction


// [s, Fs, _] = wavread("input.wav"); // _audacity.wav"); // "./signal_with_low_freq_noise_2.wav"); 

function plot_initial_things(window, signal, filtered, freq_signal, freq_filtered)

    write(%io(2), length(signal))
    write(%io(2), length(filtered))

    write(%io(2), length(freq_signal))
    write(%io(2), length(freq_filtered))

    scf(window)    
    subplot(221)
    plot2d('nl', freq_signal, abs(fft(signal_with_noise)))
    subplot(222)
    plot2d('nl', freq_filtered, abs(fft(signal_filtered)))
endfunction

load("signal_with_noise_and_filtered.sod") // Fs, signal_with_noise, signal_filtered

f1 = figure()
f2 = figure()


s = signal_with_noise(1, :);
s_len = length(s)
s_freq = (0:s_len-1)/s_len * Fs; 

s_filtered = signal_filtered(1, :);
s_filtered_len = length(s_filtered)
s_filtered_freq = (0:s_filtered_len-1)/s_filtered_len * Fs; 

plot_initial_things(f1, s, s_filtered, s_freq, s_filtered_freq)

wavwrite(signal_with_noise, Fs, 'input.wav')
wavwrite(signal_filtered, Fs, 'canonical.wav')

scf(f2)

subplot(331)
draw_freq_w_noise(s, s_freq);

cutoff =  8700/ Fs;

cutoff_low = 0 / Fs
cutoff_high = 40 / Fs

lowpass = ideal_lowpass(2048, cutoff, 0.); //  
highpass = ideal_highpass(2048, cutoff_low, cutoff_high, -1); // ideal_lowpass(256, cutoff, 0.); //  

lowpass_len = length(lowpass);
lowpass_freq = (0:lowpass_len-1)/lowpass_len * Fs; 

highpass_len = length(highpass);
highpass_freq = (0:highpass_len-1)/highpass_len * Fs; 

H_l = highpass .* lowpass
H_len = length(H_l);
H_freq = (0:H_len-1)/H_len * Fs; 


// stop();

subplot(334)
draw_ir_freq(H_l, H_freq);

ir = real(ifft(H_l));
ir_len = length(ir);

// mid = (ir_len - modulo(ir_len,2)) / 2;
shifted_ir = fftshift(ir); // cat(2, ir(mid:ir_len), ir(1:mid-1));

subplot(335)
plot2d(shifted_ir);
xlabel("Time, n", 'fontsize', 2)
ylabel("Amplitude", 'fontsize', 2)
title("Shifted IR", 'fontsize', 3)


windowed_ir = shifted_ir .* window('kr', length(shifted_ir), 16)
subplot(336)
plot2d((1:length(windowed_ir)), windowed_ir)
xlabel("Time, n", 'fontsize', 2)
ylabel("Amplitude", 'fontsize', 2)
title("Windowed Shifted IR", 'fontsize', 3)

fir_10k = windowed_ir;

conv_signal = s;

for i = 1:4
    conv_signal = conv(conv_signal, fir_10k);
end

convolved = abs(fft( conv_signal ));
conv_len = length(convolved);
conv_freq = (1:conv_len) / conv_len * Fs;


subplot(337)
plot2d('nl', s_freq, abs(fft(s)), style=color("blue"))
plot2d('nl', conv_freq, convolved, style=color("red"))

subplot(338)
part = (1:((conv_len - modulo(conv_len, 2)) / 64))
plot2d(part, conv_signal(1:length(part)))

wavwrite(conv_signal, Fs, 'conv_signal.wav');