Как вычислить коэффициенты Фурье с MATLAB

Question

Я пытаюсь вычислить коэффициенты Фурье для формы волны с использованием MATLAB. Коэффициенты могут быть вычислены с использованием следующих формул:

Т выбран равным 1, что дает омега = 2pi.

Однако у меня возникают проблемы с выполнением интегралов. Функции представляют собой треугольную волну (которая может быть сгенерирована с использованием sawtooth(t,0.5), если я не ошибаюсь), а также прямоугольную волну.

Я попытался с помощью следующего кода (для треугольной волны):

function [ a0,am,bm ] = test(numTerms) 
      b_m = zeros(1,numTerms); 
      w=2*pi; 
      for i = 1:numTerms 
        f1 = @(t) sawtooth(t,0.5).*cos(i*w*t); 
        f2 = @(t) sawtooth(t,0.5).*sin(i*w*t); 
        am(i) = 2*quad(f1,0,1); 
        bm(i) = 2*quad(f2,0,1); 
      end 
    end 

Однако это не получить в любом месте вблизи значений мне нужно. Коэффициенты b_m приведены для треугольной волны и должны быть 1/m^2 и -1/m^2, когда m нечетное чередование, начинающееся с положительного члена.

Основная проблема для меня в том, что я не совсем понимаю, как интегралы работают в MATLAB, и я не уверен, работает ли подход, который я выбрал.

Edit: Для clairify, это форма, которую я ищу, чтобы написать функцию о том, когда коэффициенты были определены:

Вот попытка с помощью FFT:

function [ a0,am,bm ] = test(numTerms) 
      T=2*pi; 
      w=1; 
      t = [0:0.1:2]; 
      f = fft(sawtooth(t,0.5)); 
      am = real(f); 
      bm = imag(f); 
      func = num2str(f(1)); 
      for i = 1:numTerms 
        func = strcat(func,'+',num2str(am(i)),'*cos(',num2str(i*w),'*t)','+',num2str(bm(i)),'*sin(',num2str(i*w),'*t)'); 
      end 
      y = inline(func); 
      plot(t,y(t)); 
    end 

Answer 1

кажется мне, что ваша проблема в том, что sawtooth возвращает mathworks documentation говорит, что:

пилообразной (т, ширина) формирует модифицированную треугольную волну, где ширина, скалярный параметр между 0 и 1, определяет точку между 0 и 2π, при котором происходит максимум. Функция увеличивается от -1 до 1 на интервале от 0 до 2π , а затем линейно уменьшается от 1 до -1 на интервале 2π шириной до 2π. Таким образом, параметр 0,5 задает стандартную треугольную волну, симметричную относительно момента времени π с амплитудой пика до пика 1. пилообразная (t, 1) эквивалентна пилообразному (t).

Так что я предполагаю, что это часть вашей проблемы.

После того, как вы ответили, я просмотрел его еще несколько. Мне кажется, что это функция quad; не очень точный! Я переделка эту проблему так:

function [ a0,am,bm ] = sotest(t, numTerms) 
     bm = zeros(1,numTerms); 
     am = zeros(1,numTerms); 
     % 2L = 1 
     L = 0.5; 
     for ii = 1:numTerms 
     am(ii) = (1/L)*quadl(@(x) aCos(x,ii,L),0,2*L); 
     bm(ii) = (1/L)*quadl(@(x) aSin(x,ii,L),0,2*L); 
     end 
     ii = 0; 
     a0 = (1/L)*trapz(t, t.*cos((ii*pi*t)/L)); 
     % now let's test it 
     y = ones(size(t))*(a0/2); 
     for ii=1:numTerms 
     y = y + am(ii)*cos(ii*2*pi*t); 
     y = y + bm(ii)*sin(ii*2*pi*t); 
     end 
     figure; plot(t, y); 
    end 

    function a = aCos(t,n,L) 
     a = t.*cos((n*pi*t)/L); 
    end 

    function b = aSin(t,n,L) 
     b = t.*sin((n*pi*t)/L); 
    end 

И тогда я назвал это нравится:

[ a0,am,bm ] = sotest(t, 100); 

и я получил:

Сладость !!!

Все, что я действительно изменился был от quad к quadl. я понял, что из с помощью trapz, которые хорошо работали, пока вектор времени я использовал не имел достаточного разрешения, что привело меня не поверить, что это числовая вопрос, а не что-то фундаментальна. Надеюсь это поможет!

Answer 2

Чтобы устранить неполадку в коде, я бы построил функции, которые вы используете, и исследуйте, как их выполняет выборка. Возможно, вы несете их сэмплирование, поэтому убедитесь, что минимальный размер шага меньше, чем период функции, по меньшей мере, на коэффициент 10.

Я бы предложил использовать БПФ, которые встроены в Matlab. Мало того, что БПФ является наиболее эффективным методом вычисления спектра (это n * log (n) зависит от длины n массива, тогда как интеграл от n^2 зависит), он также автоматически даст вам частотные точки которые поддерживаются вашими (равноотстоящими) данными времени. Если вы сами вычислите интеграл (может потребоваться, если точки данных не равномерно распределены), вы можете рассчитать частотные данные, которые не разрешены (более близкое расстояние, чем 1/по интервалу во времени, то есть за пределами «предела Фурье»).