Fir_filter-matlab-code.md

@@ -0,0 +1,146 @@
+function custom_filter_design_tool_4()
+    f = input('Enter the frequency (in Hz): ');
+
+    if f >= 1000
+        [Fpass, Fstop, Wpass, Wstop, N, dens] = set_values_forxxx(f);
+        filter_type = 'firpm';
+    else
+        [Fpass, Fstop, Wpass, Wstop, N, dens] = set_values_yy(f);
+        filter_type = 'firls';
+    end
+
+    Fs = 40000;
+
+    % Filter design based on parameters
+    if f < 1000
+        b = firls(N, [0 Fpass Fstop Fs/2]/(Fs/2), [1 1 0 0], [Wpass Wstop]);
+    else
+        b = firpm(N, [0 Fpass Fstop Fs/2]/(Fs/2), [1 1 0 0], [Wpass Wstop], dens);
+    end
+
+    % Display the filter type
+    disp(['Filter Type: ' filter_type]);
+
+    % Initialize the frequency step
+    step_size = input('Enter the frequency step size (50 or 100 Hz): ');
+    if step_size ~= 50 && step_size ~= 100
+        error('Invalid step size. Please enter either 50 or 100 Hz.');
+    end
+
+    % Compute magnitude at the given frequency and iterate until the nearest multiple of the step size
+    while mod(f, step_size) ~= 0
+        if f > (f - mod(f, step_size))
+            f = f - step_size;
+        else
+            f = f + step_size;
+        end
+    end
+
+    while true
+        [magnitude_dB,last_Fpass,Fstop_last, Wpass_Last,Wstop_last, c] = compute_magnitude_at_frequency(b, f, Fs, N, Fpass, Fstop, Wpass, Wstop, dens);
+        disp(['Magnitude at ' num2str(f) ' Hz: ' num2str(magnitude_dB) ' dB']);
+        disp(['Last value of Wpass: ', num2str(Wpass_Last)]);
+        disp(['Last value of Wstop: ', num2str(Wstop_last)]);
+        disp(['Last value of Fstop(Hz): ', num2str(Fstop_last)]);
+        disp(['Last value of Fpass(Hz): ', num2str(last_Fpass)]);
+
+        % Plot magnitude response
+        figure;
+        freqz(c, 1, 1024, Fs);
+        filename = sprintf('filter_coefficients_6_%dHz_%s.txt', f, filter_type);
+        save_filter_coefficients(filename, c);
+
+        % Update frequency for the next iteration
+        f = f - step_size;
+        if f <= 0
+            break;
+        end
+    end
+end
+
+function save_filter_coefficients(filename, c)
+    fid = fopen(filename, 'wt');
+    fprintf(fid, '%f\n', c);
+    fclose(fid);
+    disp(['Filter coefficients saved to ' filename]);
+end
+
+function [Fpass, Fstop, Wpass, Wstop, N, dens] = set_values_forxxx(f)
+    % For frequencies above 1 kHz
+    Fpass = f - 100;
+    Fstop = f + 100;
+    Wstop = 755;
+    Wpass = 2;
+    N = 8;
+    dens = 20;
+end
+
+function [Fpass, Fstop, Wpass, Wstop, N, dens] = set_values_yy(f)
+    % For frequencies below 1 kHz
+    Fstop = f + 10;
+    Fpass = f - 10;
+    Wstop = 1;
+    Wpass = 1;
+    N = 22;
+    dens = 20;
+end
+
+function [magnitude_dB, last_Fpass, Fstop_last, Wpass_Last,Wstop_last, c] = compute_magnitude_at_frequency(b, f, Fs, N, Fpass, Fstop, Wpass, Wstop, dens)
+    while true
+        if f < 1000
+            c = firls(N, [0 Fpass Fstop Fs/2]/(Fs/2), [1 1 0 0], [Wpass Wstop]);
+        else
+            c = firpm(N, [0 Fpass Fstop Fs/2]/(Fs/2), [1 1 0 0], [Wpass Wstop], dens);
+        end
+
+        [H, freq] = freqz(c, 1, 1024, Fs);
+        [~, idx] = min(abs(freq - f));
+        magnitude_dB = 20*log10(abs(H(idx)));
+
+        % Adjust Fpass and Fstop based on the magnitude
+        if f < 1000
+            if magnitude_dB < -48
+                Wpass = Wpass + 1;
+            elseif magnitude_dB > -45
+                Fstop = Fstop - 10;
+                Fpass = Fpass - 10;
+                Wstop = Wstop + 10;
+            else
+                break; % Exit loop if magnitude is within desired range
+            end
+        else
+            if magnitude_dB < -46 && f >= 3000
+                Fpass = Fpass + 100;
+                Fstop = Fstop + 100;
+                Wstop = Wstop - 10;
+            elseif magnitude_dB > -40 && f >= 3000
+                Fpass = Fpass - 100;
+                Fstop = Fstop - 100;
+                Wstop = Wstop + 10;
+            elseif magnitude_dB < -47 && f < 3000
+                Fpass = Fpass + 10;
+                Fstop = Fstop + 10;
+                Wstop = Wstop - 1;
+            elseif magnitude_dB > -40 && f < 3000
+                Fpass = Fpass - 10;
+                Fstop = Fstop - 10;
+                Wstop = Wstop + 10;
+            else
+                break; % Exit loop if magnitude is within desired range
+            end
+        end
+    end
+
+    % Return the last value of Wpass
+    % last_Wpass = Wpass;
+    % Fstop_last = Fstop;
+    % Fpass_Last = Fpass;
+    % Wstop_last = Wstop;
+        last_Fpass = Fpass;
+        Fstop_last = Fstop;
+        Wpass_Last = Wpass;
+        Wstop_last = Wstop;
+    Hd = dsp.FIRFilter('Numerator', c);
+    
+    freqz(Hd);
+end