texts:
I am new to DSP, and currently trying to implement a simple polyphase decimation program in matlab.

M is decimation factor, x is input (for now it is a square wave 8096 samples), y_polyDec is result from polyphase, output is result from regular filter then decimation.

code :

FilterBank_LPF = fir1(127, 1/M);
y_polyDec = zeros(8096/M, 1);

for i = 1:M

    p(:, i) = FilterBank_LPF(i:M:end);
    x(:, i) = input(M-i+1:M:end);
    y_polyDec = filter(p(:, i), 1, x(:, i))+ y_polyDec;

end

out_LPF = filter(FilterBank_LPF, 1, input);
output = downsample(out_LPF, M);

%===========================================

If I compare "output" and "y_polyDec" (either comparing the decimated sample values or frequency spectrum using fvtool) I notice there are differences in the values. I was wondering if my code is wrong, or am I still missing something.

Thank you


texts:
Why is the Fourier transform valid only for absolutely integrable signals?

For example, why can't we do the Fourier transform of exponential order functions?


texts: Effect of down mixing on complex demodulation method

texts:
As I know, IQ demodulation (complex demodulation) consists of three steps:

Down mixing
X(complex demodulation)(t)=X(real signal)(t).{exp[(-i2π(f_demodulation).t)]}

low pass filtering

decimation

would you please explain what consequence would effect the demodulation procedure by implementing down mixing conversion in step1?
As a second question, when do we use complex demodulation instead of regular methods for demodulating the signal to extract specific frequencies?
Regards
Mahsa


texts: Relationship between the autocorrelations of X(t) and X(nt)

texts:
Defining:

$X(t)$ WSS random process with autocorrelation function $R_{X}(\tau) = \mathbb{E}[X(t)X(t+\tau)]$.

$Y[n] = X(nT)$ (sampling of $X$ at a rate $\frac1T$) with autocorrelation function $R_Y(\tau) = \mathbb{E}[Y[n]Y[n+\tau]]$.

Which is the relationship between $R_X(\tau)$ and $R_Y(\tau)$?


texts: Has anyone used Siglib? Answer a question on the FFTr function?

texts:
I'm hoping someone - anyone can answer a question on the Siglib library call to perform an FFT on data.

There's no examples anywhere I can find that provide a real application using data from an ADC, so that's where my wheels are spinning.
The examples from Numerix have only canned data, generated by the function SDA_SignalGenerate. Specifically, in the FFT example, SDA_SignalGenerate is passed a value of $0.15625$ for frequency. Can anyone provide some insight as to what this represents? The FFT_LENGTH of the FFT function is $512$, but I do not see what the SDA_SignalGenerate frequency of $0.15625$ means, or what it is normalized to.

Thanks for any help anyone can provide.


texts: How can understand periodicity of a Signal from frequency domain representation?

texts: Is it possible to say a signal is periodic from its frequency domain representation? A periodic signal is sum of its sinus and cosinus. Frequency translation of sinus and cosinus functions are symmetrical impulses. Is it correct to say "a signal is periodic if its frequency domain representation includes symmetrical inpulses"?


texts: What is the relation between the domain of Fractional Fourier transform and time and frequency?

texts:
This refers to the link where it clearly writes that:

  The FrFT argument $u$ is neither a spatial one $x$ nor a frequency
  $\xi$. We will see why it can be interpreted as linear combination of
  both coordinates $(x,\xi)$

But I did not see the relation between the domain of Fractional Fourier transform ($u$) and space and frequency $(x,\xi)$. I can understand that domain '$u$' is a linear combination of '$(x,\xi)$' but my question what is the relation? Can anyone point out the exact relation?


texts: What is the correct way to do Short term Fourier transform and extract the phase information from local sections of a signal?

texts:
I took a triangular voltage signal and taking it's time derivative I obtained the current through a pure capacitor.
Then I took different portion of voltage and current signals, took their Fourier transform (as if I am taking the STFT with rectangular gate function) and obtained the magnitude & phase from FFT. I obtained the sampling frequency by the formula- Fs = N/(t_final-t_initial) where t_initial & t_final are the time of the starting point and ending point of the portion of signal I choose.

For pure capacitor, I expect $+90^0$ phase shift of each current components with respect to corresponding voltage components. But depending on the portion of signal I choose, sometimes I find $+90^0$, sometimes $+270^0=-90^0$, sometimes $+60^0$!!! In some cases I don't see any meaningful result (see 3rd figure from left)! In addition, different frequency component shows different phase relation.

So, am I doing any mistake? How to obtain correct phase information at every portion of the signal?

MATLAB code:

clc;
clear;

%% Generate the signal to be processed
C=10e-6;       %10 micro-farad

t_full=[1:500]*1e-6;
v1=1*t_full(1:end/2);  %1st Ramp of the triangle
v2=v1(end)-1*(t_full(end/2+1:end)-t_full(end/2)); %2nd Ramp of the triangle
v_full=[v1 v2];  %Triangular voltage signal

dv=gradient(v_full); dv_dt=dv./(t_full(2)-t_full(1));
i_full=C*dv_dt;  %Capacitor current =C*dv/dt

%% Choose the portion of the full signal for STFT
%Total signal has 500 points.

st=40;  %starting point of the signal portion
np=125; %ending point of the signal portion

t=t_full(st:st+np);
v=v_full(st:st+np);
i=i_full(st:st+np);

%% Calculate sampling frequency
N = length(t);                % Length of signal under study
t_final=t(N);t_initial=t(1);
Fs = N/(t_final-t_initial);   % Sampling frequency

%% Obtain the FFT
V = fft(v)/N;
I = fft(i)/N;

mag_V_full = abs(V);
mag_V = mag_V_full(1:N/2+1);mag_V(2:end-1) = 2*mag_V(2:end-1); %FFT is symmetric, taking positive frequency only
phi_V_full=angle(V);phi_V = phi_V_full(1:N/2+1);

mag_I_full = abs(I);
mag_I = mag_I_full(1:N/2+1);mag_I(2:end-1) = 2*mag_I(2:end-1); %FFT is symmetric, taking positive frequency only
phi_I_full=angle(I);phi_I = phi_I_full(1:N/2+1);

f = Fs/N*(0:(N/2));

del_phi=(phi_V-phi_I)/pi*180;  %Phase shift between voltage & current component

%% steming/ploting result

subplot(4,1,1)
[hAx,hLine1,hLine2] =plotyy(t,i,t,v);
xlabel('time (s)')
ylabel(hAx(1),'Current (A)'); % left y-axis
ylabel(hAx(2),'Voltage (V)'); % right y-axis

subplot(4,1,2)
stem(f/1e3,mag_V)
ylabel '|V|'
xlim([-0.05 50])

subplot(4,1,3)
stem(f/1e3,mag_I)
xlim([-0.05 50])
ylabel '|I|'

subplot(4,1,4)
stem(f/1e3,del_phi)    %We need to see the phase shift of the dominating frequency component.
xlim([-0.05 50])
ylabel 'angle of Z=(Phi_V-Phi_I)'
xlabel 'Frequency (kHz)'


texts: UPDATE : How to continue Computing the Convolution

texts:
$$
x(t):=\begin{cases}
1&\text{if $0<t<T$}\\ \\
0&\text{if otherwise}
\end{cases}
\qquad\text{and}\qquad
h(t):=\begin{cases}
t&\text{if $0<t<2T$}\\ \\
0&\text{if otherwise}
\end{cases}
$$
Compute $(x*h)(t)$.
My approach:
We first note that :$$\displaystyle x(\tau)=\begin{cases}
1&\text{if $0<\tau<T$}\\ \\
0&\text{if otherwise}
\end{cases}\qquad\text{and}\qquad\displaystyle h(t-\tau)=\begin{cases}
t-\tau&\text{if $0<t-\tau<2T$}\\ \\
0&\text{if otherwise}
\end{cases}$$ We observe that :
$$
0<\tau<T
$$
$$
t-2T<\tau<t
$$
Thus, $(x*h)(t)\in\mathcal{R}[\max(0,t-2T),\min(t,T)]$ (i.e. Riemann integrable). Therefore :
\begin{align*}
(x*h)(t)&:=\int_{-\infty}^{\infty}x(\tau)h(t-\tau)\;\text{d}\tau\\
&=\int_{\max(0,t-2T)}^{\min(t,T)}t-\tau\;\text{d}\tau\\
&=\left. t\tau-\frac{\tau^{2}}{2}\right|_{\max(0,t-2T)}^{\min(t,T)}\\
&=\left[t\min(t,T)-\frac{\min^{2}(t,T)}{2}\right]-\left[t\max(0,t-2T)-\frac{\max^{2}(0,t-2T)}{2}\right]\\
&=t[\min(t,T)-\max(0,t-2T)]-\frac{1}{2}[\text{min}^{2}(t,T)-\text{max}^{2}(0,t-2T)]
\end{align*}

Problem: I want to further solve this to obtain the final solution as a piecewise function (if it has such form)

I present to you the MATLAB simulation of the convolution

First one corresponds to $T=1$ and the second one corresponds to $T=2$. Here is the algorithm used:
t = linspace(0,10,1000);   % assume a time span for "t"
T=1;
m = min(t,T);              % find min of (t-1 , 5) for each "t"
M = max(0,t-2*T);            % find max of (t-3 , 3) for each "t"
vind = M < m;              % find "valid" range of "t"

y = t.*(m-M)-0.5*(m.^2-M.^2);% EVALUATE the CONVOLUTION
y = y.*vind;               % force non-valid range to zero.

figure,plot(t,y);         % DISPLAY:
title('convolution is :');

Credits to $\text{@Fat32}$ for the code in his answer in the link provided in the comment


texts: Coefficients of the Fourier Series expansion of physical signals having finite power

texts:
Question 2.7. Show that for all periodic physical signals that have finite power, the coefficients of the Fourier series expansion $
x_n $ tends to 0 as $ n \to \infty $

I have computed $ |x_n|^2 $ by multiplying $ x_n $ with $ x_n^*$. However, as per my calculations, $ |x_n|^2 = P_x $ i.e. power of $
x(t) $, which is a constant and independent of $n$. I know that this is fundamentally incorrect because $ |x_n| $ should vary with $n$. If anyone can guide me about where I went wrong that would be very helpful. Thanks and regards.

  MY SOLUTION:

Please refer to the solution from the half of the first page.


texts: Normalizing a vector with Mel Cepstrum Coefficients, Delta Coefficients and Log Energies for both

texts:
This is my very first signal processing exercise so I'm very new to this space. That being said the goal is to do isolate voice recognition in identifying someone saying yes and no.As the title suggests I have a 1x26 vector that looks as follows:

  Columns 1 through 10

  (P)160.2504   66.7557  -35.9298  144.1598  -11.9133  -16.2184  -26.4730   78.0673   10.2075  -39.4318

  Columns 11 through 20

  -46.6702    5.9283  -42.9540    (P)0.2614    0.0382   -1.3315    0.1326   -1.4256    0.1936    0.4081

  Columns 21 through 26

   -1.4431   -1.6916   -1.1668   -0.2901   -0.4499   -0.4309

Where (P) denotes the log energy for the Mel Cepstrum Coefficients and Delta Coefficients respectively. This being said I would like to know how to normalize these values but how? I'm not sure what exactly to divide the energies by, let alone the coefficients. Like I said first project so a really go post guiding me would be REALLY appreciated!


texts:
This question is about SNR in the context of EEG. (a related question  about SNR)

I am interested in calculating The SNR of an ERP. My motivation is: To calculate the "signal-to-noise ratio" and from this to extract the information (how?) on how many trials are needed to construct a stable ERP.

my Main question is: how do I calculate and SNR for such scenario  ?

some more explation about EEG and ERP are enclosed:
EEG:

  The on-going electrical activity of the brain measured from scalp electrodes is called the electroencephalogram or EEG.  within the EEG is a signal which is more revealing about information processing in the brain. This signal can be obtained by time-locking the recording of the EEG to the onset of events such as a person reading a word on a computer screen, listening to a musical note played on an instrument, or viewing a picture in a magazine. The resulting activity is called an "event-related potential" (ERP)

SNR:

  The number of trials necessary to obtain an ERP depends on a number of factors, the most important being the "signal-to-noise ratio", that is, the relative size of the signal (the ERP) relative to the size of the noise (the background EEG). In cognitive experiments 30 to 50 stimulus presentations are typically required to obtain a good clean average ERP.

from: http://neurocog.psy.tufts.edu/images/erp.htm


texts: Where to get transtion matrix for Kalman filter?

texts:
I am trying to understand the Kalman filter, and I am struggling to find out how to choose the transition matrix (denoted as $\mathbf F(k)$ or $\mathbf A(k)$. This matrix is used to update the weight error covariance matrix as follows

$$\mathbf P(k+1) = \mathbf F(k+1) \mathbf P(k) \mathbf F^\mathrm T (k) + \mathbf Q$$

In the paper On the Intrinsic Relationship Between the Least Mean Square and Kalman Filters it is description of "The Kalman filter for deterministic states.", what does not use the part with prediction, so it is without the transtion matrix.

But there is also the general Kalman filter, with the prediction according to transtion matrix. I have not understand how this matrix should be obtained.

I quess, that this transition matrix is from state space formulation of the model. But what should I do if I do not have a model? Is this mean that the general Kalman filter is usable only If I have exact model of the system?


texts: What is the "bilinear interpolation kernel" in personlab paper?

texts:
please excuse my ignorance in computer graphics, but what is this bilinear interpolation kernel in the personlab paper page 6 equation (1)?

Here it is:
$$
h_k(x) = \frac{1}{\pi R^2}\sum_{i=1:N}p_k(x_i)B(x_i+S_k(x_i)-x)
$$
where:

$x$, $x_i$ are image coordinates (2-D vectors).
$h_k \colon \mathbb{R}^2 \mapsto \mathbb{R}$ localized heatmaps obtained by a sort of Hough voting.
$S_k \colon \mathbb{R}^2 \mapsto \mathbb{R}$ offest vectors
$p_k \colon \mathbb{R}^2 \mapsto \mathbb{R}^2$ heatmaps
$B(.)$ denotes the bilinear intepolation kernel.

  This is a
  form of Hough voting: each point i in the image crop grid
  casts a vote with its estimate for the position of every keypoint, with the vote being weighted by the probability that
  it is in the disk of influence of the corresponding keypoint

I guess bilinear interpolation kernel is a class of function and not an actual one (I may be wrong).

But if it is bilinear why does it only take one argument in the equation?

Maybe this is bilinear 1D interpolation, intepolating between the coordinates of 2-D vector?

Otherwise this formula doesn't make sense...


texts: Understanding noise removal method using wavelets

texts:
I am trying to understand how wavelet transform can be used to denoise a time series or signal and how to plot the scalogram image. My signal has a lot of fluctuations and as such I am finding it difficult to denoise. Morevoer, to plot the scalogram I need to know the frequency. I don't know what is the frequency for this particular kind of time series obtained from a dynamical system of the form: Logistic Map given by:
$$x[n] = 4\big(x[n-1]\big)\big(1-x[n-1]\big)$$
Systems similar to this type of dynamical systems are the Lorenz, Mackey-Glass. Can somebody please help:

How to properly denoise the signal? As observed, from the plot the denoised signal denoised does not look exactly the same as the clean signal x (black dotted line), so what are other parameters or wavelet types that I could use and how to decide which ones to use. Is there a rule of thumb?

What is the sampling and nyquist frequency for this kind of signal and

how to plot the scalogram image: I used wt() to obtain the wavelet coefficients. After that how to plot the image of scalogram so that X axis is time and Y axis is Frequency?
 x(1) = 0.1; % initial condition (can be anything from 0 to 1)
 M = 50; %number of data points (length of the time series)
 for n = 2:M, % iterate
     x(n) = 4*x(n-1)*(1-x(n-1));
 end

 %add noise
 x_noise = awgn(x,10,'measured');

 %denoise using wavelet
 denoised = wdenoise(x_noise, 3,'Wavelet','db3',...
     'DenoisingMethod','Bayes',...
     'ThresholdRule','Median',...
     'NoiseEstimate','LevelIndependent');
 figure
 plot(x_noise)
 axis tight
 hold on
 plot(denoised,'r')

 fb = cwtfilterbank('SignalLength',M);
 [cfs,frq] = wt(fb,denoised);


texts:
In many applications such as estimation theory, when we need to estimate a parameter then we usually consider in presence of white gaussian noise of zero mean and some standard deviation. During Maximum likelihood estimation, we also use this assumption. So, my question is -

Do we consider noise to be uncorrelated or correlated in estimation?
What is the difference between correlated and uncorrelated noise and its significance
Why do we consider such a property of correlated or uncorrelated in estimation and when we say that measurement noise is gaussian.


texts:
I'm developing a PID controller to control the position of a laser beam. The position is measured using a Position Sensitive Device (PSD). I modeled the plant according to the manufacturer specs and I'm using Matlab/Simulink with PID tuner and embedded coder to generate the code for an F28335 micro-controller.

As long as I use a PI controller, the results are consistent with the model, but as soon as I introduce a differentiator and try to obtain a sharper step response, the system becomes unstable, with a lot of oscillation and saturation.

My guess is that the problem is the noise in the feedback sensor, so I tried to filter both the feedback signal and the error to no avail. Is there any other reason why the real system becomes so unstable? Any suggestions?

Edit:

This is the compensator formula provided by the Simulink PID block:


texts: Downsampling - how many extra samples to include

texts:
I want to downsample a finite section of an equidistant time-series
$$x_i, \qquad i=0\dots n-1$$
 by doing an FFT and dropping the upper part of the spectrum. Since the FFT looks at $x_i$ as a cyclic series, it will show edge effects (ringing) near $i=0$ and $i=n-1$ caused by for example a jump between $x_{n-1}$ and $x_0$ in the cyclic view. Hence I need to include some extra samples left of $x_0$ and right of $x_{n-1}$, i.e. take the indices
$$i=-m\dots n-1+m$$
into the FFT, where the corresponding samples $x_{i<0}$ and $x_{i>n-1}$ are either accessible values (measurements) or respresent a reasonable extrapolation (e.g. zero-padding, $x=0$). The extra samples are dropped again once the downsampling is completed.

Now the question is, how many samples should I add? It is clear that the 'sinc' function, which represents the antialiasing cutoff (and hence, the edge ringing) is infinitely extended. So, I can, of course, not expect to obtain an exact solution from adding a finite set of samples. On the other hand, some quick examples show that adding just one sample on either side is not enough in most cases, since this only goes one half swing of the 'sinc' function towards the boundary values.

What heuristics could lead me to choose a reasonably appropriate number of extra samples (probably by specifying an accuracy parameter)?

The image above shows the result of a downsampled time-series that originally started at 1, but has been extended by a bunch of zeroes to the left. Since the ringing is not only determined by the sample rate but also by the height of the original step (in this case =1), it is not so easy to determine when the ringing (sinc) has dropped 'enough'. Imagine the step scaled by 1000, then the ringing is also scaled by 1000.


texts: Is the behaviour of the following filter plausible?

texts:
I've designed a filter based on this equation

$$
E(u) = \frac{1}{2}\int_{\Omega} \left[(u-f)^2 + g(\lVert \nabla u \rVert)\right]dxdy
$$

Euler-Lagrange equations in this case are (should be) given by

$$
u - f - \frac{1}{2}\nabla^T \cdot \left( \frac{\nabla u}{\lVert \nabla u\rVert}g'(\lVert \nabla u\rVert)\right) = 0
$$

For given $\lambda > 0$ if I chose $g(x) = e^{-x/\lambda}$ the EL equations are
$$u - f + \frac{1}{2\lambda}\nabla^T \cdot \left( \frac{\nabla u}{\lVert \nabla u\rVert}e^{-\lVert \nabla u\rVert/\lambda}\right) = 0$$

Continuous gradient descent in this case is given by

$$
u_t = f - u - \frac{1}{2\lambda}\nabla^T \cdot \left( \frac{\nabla u}{\lVert \nabla u\rVert}e^{-\lVert \nabla u\rVert/\lambda}\right)
$$

Basically my idea would be... if the gradient magnitude is very high then $u = f$ (i.e. I'd like to preserve strong edges) otherwise (small variations) $u$ should be tending to some constant, or anyway smaller variations.

Now... This is my implementation in opencv

void g(InputArray input, OutputArray output) {
    //cv::exp(-im,out);
    float lambda = 0.0000001;
    Mat rval, arg = input.getMat();
    //cv::exp(-arg / lambda2, tmp);
    cv::exp(-arg/lambda,rval);
    output.assign(rval);
}

void dg(InputArray input, OutputArray output) {
    Mat rval, tmp, arg = input.getMat();
    g(arg,tmp);
    float lambda = 0.0000001;
    rval = (-1.0/lambda)*tmp;
    output.assign(rval);
}

void evolveImage(const cv::Mat & input, cv::Mat & output, void (*g)(cv::InputArray, cv::OutputArray), void(*dg)(cv::InputArray, cv::OutputArray), const int & k_max) {
    int X = input.cols;
    int Y = input.rows;

    Mat imfloat;
    input.convertTo(imfloat, CV_32FC1);
    Mat u;
    imfloat.copyTo(u);
    cv::Mat & f = imfloat;

    double minval, maxval;

    minMaxLoc(u, &minval, &maxval);
    Mat gradientMagnitude;

    Mat dx, dy;
    Mat gprime, a, b;
    for (auto k = 0; k < k_max; ++k) {

        Sobel(u, dx, CV_32F, 1, 0, 3); //Estimation differences in x
        Sobel(u, dy, CV_32F, 0, 1, 3); //Estimation differences in y
        cv::sqrt(dx.mul(dx) + dy.mul(dy),gradientMagnitude);
        dg(gradientMagnitude, gprime);
        minMaxLoc(dx, &minval, &maxval);
        minMaxLoc(dy, &minval, &maxval);
        a = gprime.mul(dx);
        b = gprime.mul(dy);
        Sobel(a, dx, CV_32F, 1, 0, 3);
        Sobel(b, dy, CV_32F, 0, 1, 3);

        u = f + 0.5f*(dx + dy);

    }

    minMaxLoc(u, &minval, &maxval);
    u.convertTo(output, CV_8UC1, 255 / (maxval - minval), -minval * 255 / (maxval - minval));
}

I have this messy function that can be called by a main if you do cut and paste

int continuousGradientDescentExample(int argc, char **argv) {
    if (argc != 2) {
        std::cout << "Need filename..." << std::endl;
        return 1;
    }

    Mat im = imread(argv[1], CV_LOAD_IMAGE_GRAYSCALE);

    namedWindow("Test", WINDOW_AUTOSIZE);
    imshow("Test", im);
    waitKey(0);

    Mat imfloat;
    im.convertTo(imfloat, CV_32FC1);
    Mat u(imfloat.size(), imfloat.type());

    int & X = im.cols;
    int & Y = im.rows;

    evolveImage(imfloat, u, g, dg, 50);

    std::cout << type2str(u.type()) << std::endl;

    Mat output(Size(2 * X, Y), CV_8UC1);
    im.copyTo(output(Rect(0, 0, X, Y)));
    u.convertTo(u, CV_8UC1);
    u.copyTo(output(Rect(X, 0, X, Y)));

    namedWindow("Original vs Anisotropic diffusion", WINDOW_AUTOSIZE);
    imshow("Original vs Anisotropic diffusion", output);

    waitKey(0);

    return 1;
}

The code is working (i.e. doesn't crash, you can try it if you want). But the thing is no matter how I choose $\lambda$ in the evolveImage the matrices $a,b$ are always $0$ which implies the update of $u$ doesn't produce any effect. To me the equation makes sense, but I can't understand why there's no difference at all.

Can anyone give me an insight?


texts:
[
Here are two figures from the same signal, the 2nd one(B) is after digital predistortion at transmitter end. How can you describe it? What are the improvements and advantages of this technique?
Much appreciate your answers.


texts:
I am searching for global purpose Dataset to test some image registration algorithm.. May you please give advice?
I am looking for color or gray-scale dataset. I prefer it not to be medical or binirized.
thanks

EDIT:

I found this dataset: http://www.robots.ox.ac.uk/~vgg/data/data-aff.html
it seems good but it contains just few images (6 in each set)


texts: Compare two spoken words with MFCC and DTW

texts:
Processing of the speech signal is done using Matlab. I want to use MFCC (Mel-frequency cepstral coefficients) is implemented using HTK and DTW (dynamic time warping v2.1) algorithms to make a speech comparison script. I calculated MFCC and got 13 coefficients for each frame.

How do I calculate DWT from Applying DTW on MFCC frames?
DTW gives me a distance between two patterns. How do I determine if they match or not?

How do I edit this code?
this is the code I am using:
function d=dtw(s,t,w)

% s: signal 1, size is ns*k, row for time, colume for channel
% t: signal 2, size is nt*k, row for time, colume for channel
% w: window parameter
% if s(i) is matched with t(j) then |i-j|<=w

if nargin<3
    w=Inf;
end

ns=size(s,1);
nt=size(t,1);

if size(s,2)~=size(t,2)
    error('Error in dtw(): the dimensions of the two input signals do not match.');
end

w=max(w, abs(ns-nt)); % adapt window size

%% initialization

D=zeros(ns+1,nt+1)+Inf; % cache matrix
D(1,1)=0;

%% begin dynamic programming

for i=1:ns

for j=max(i-w,1):min(i+w,nt)

    oost=norm(s(i,:)-t(j,:));

    D(i+1,j+1)=oost+min( [D(i,j+1), D(i+1,j), D(i,j)] );

end

end

d=D(ns+1,nt+1);


texts: Finding the maximum phase and the maximum frequency deviation of an angle-modulated signal

texts:
I have an angle-modulated signal:
$$
s(t)=10\left[\cos\left(10^8\pi · t + 5\sin\left(2\pi\cdot10^3·t\right)\right)\right]
$$
How would I go about finding the maximum phase and the maximum frequency deviation?


texts: Dynamic range compression - Need pseudo algorithm for normalizing a signal

texts:
I have a signal that is NOT an audio signal (but looks like an audio signal - probably you can even "play" it). I would like to amplify the signal in regions where its amplitude is low. In very rare occasions the signal amplitude may be too high so I want to decrease it a bit. At the beginning and at the end of the 'track' I have some noise.

I think what what I want is called "Dynamic range compression". Is this correct? How I do this programmatically? I need some pseudo-code.


texts: Difficulties in understanding mutual information concept

texts:
What does mutual information (MI) convey? If 2 signals are independent then MI is zero; What does this imply and mean in the case of mutual information of error, and in general definition of MI and entropy.
What is the interpretation if MI of error is decreasing or getting minimized?


texts:
I am performing an analysis of an output signal in the interval $t<0$ from the operational amplifier given in the following image

where $v_{C_1}(0^-)$ and $v_{C_2}(0^-)$ are the initial conditions to be computed of two capacitors. As you can see, the solution is just below the image.

What I can't understand is why all the measured voltages are $1V$ without knowing the values of the resistances and just knowing that $V_-=V_+$. For sure I am missing something that I don't know. Can someone give me a hint? Thank you.


texts: Proper generation of Complex white Gaussian noise using Matlab

texts: I am not comfortable in using the Matlab's awgn() function as it is hard to understand what is actually going on. So, I wanted to confirm if the following way is correct to generate noisy samples of a particular snr of 30 dB and a particular variance for data in complex domain.


texts: Is ergodic in mean a property defined only for WSS stochastic processes?

texts:
I understand the definition of a random process $X(t)$ being ergodic in mean (first-order ergodic) is that the expectation of the sample mean $<u_X>_T=\frac{1}{T} \int_{\frac{-T}{2}}^{\frac{T}{2}} X_w(t)dt $ should converge to the ensemble mean as $T \rightarrow \infty$. Can we define ergodic in mean property for processes which do not have constant mean.

I know that the property of being ergodic in ACF, second-order ergodic, requires WSS assumption.

My question is motivated by the fact that the requirement of WSS is not given in the wiki article describing ergodic process but used in the example illustrating the same.


texts: How to apply an offset to a complex exponential without chopping it off?

texts:
I am trying to encode a bitstream as offset in a complex exponential, but when I do, the sinusoids get chopped off. Do I need to apply a lowpass to the bitstream or is there a more elegant way?
bits_t bit;
double t, cycle;
double pi = 4 * atan(1.0);
doubl x, y;
double complex z;
short ch[2];
int f = CARRIERFREQ;
int fs = SAMPLERATE;

z = cexp(-I * 2 * pi * f * n/fs);
z = (z + (bit+I*bit))/2;
x = creal(z);
y = cimag(z);

The generation of the bistream is not included in the code. Just take the variable bit given as the current bit.

I basically want them to look like this


texts: Is there any practical application for performing a double Fourier transform? ...or an inverse Fourier transform on a time-domain input?

texts:
In mathematics you can take the double derivative, or double integral of a function. There are many cases where performing a double derivative models a practical real-world situation, like finding the acceleration of an object.

Since the Fourier transform takes a real or complex signal as an input, and produces a complex signal as an output, there is nothing stopping you from taking that output and applying the Fourier transform a second time... Are there any practical uses for doing this? Does it help to model some complex real-world situations?

With the same logic, nothing would stop you from taking the inverse Fourier transform of your original time-domain input signal... would this ever be useful? Why or why not?


texts: Removing mean from signal massively distorts FFT

texts:
I have a series of signals sampled at 1200 Hz that are 200-300 samples long. I'd like to analyze them in the frequency domain in 1 Hz bins in MATLAB. I'm primarily interested in the frequencies below about 20 Hz, since this is where virtually all of the signal power is contained.

However, I'm finding conflicting opinions on whether it's appropriate to subtract the mean value from the signal before I perform the FFT function.  When I don't subtract the mean, I get a strong 0Hz component, which makes sense given the shape of the original signal. But some sources say you should zero-mean the data so you don't get distortion of the low frequencies, which are the ones I care about.  If I do subtract the mean, the FFT is markedly different, and appears to contain massive sinc-like distortions.

Here is the original signal in blue, and the zero-mean version in red, created using detrend(mySignal, 'constant').

As expected, the detrended signal has a mean of zero. But when I do the FFT and look at the spectral distribution, this is what I get:

Note that X-axis is frequency and I've truncated the axes to the first 50 frequency bins (0-49 Hz).

Is it inappropriate to zero-mean a signal like this? And what's the meaning of this distortion in the zero-meaned signal, if anything? Any help is much appreciated!


texts:
I want to build a detector of sorts. Say I have a bunch of signals and they all share some patterns, like some peaks in frequency or something more complicated that I can't bother to calculate by hand. I do however have plenty of these signals, so maybe I could "model" statistically, and more importantly, numerically through Matlab, their behavior. I'm familiar with the concepts of auto correlation and spectral density, but I'm not sure how to infer them from of the samples. I know Matlab has a function called autocorr that calculates the autocorrelation of the sample, but I'm not sure I understand how it works. In any case, I don't know how to get the best autocorrelation model out of all the autocorrelations calculated for each sample. Taking an average seems a bit too easy and I'm not sure how it would work.

Thank you!


texts: Is there any example or existing resource for using matlab in audio comparison (cross correlation)?

texts: I am new to MATLAB. I have the audio files A and B, B is recorded using a mic while playing audio file A. The length of file B is different from A. A does not always start from 0 second. My program is supposed to identify, file B from file A (giving a set of files A, X, Y, Z, etc. the program can figure out that the most correlation is A) Is there any existing resource or example for my problem? Or any suggestions for me to get start? Appropriate for your any comments and thanks!


texts: How to deconvolve matrix when a model of noise exists?

texts:
I have a matrix of N rows of time-series data. There is a specific noise contaminating measurement of the data that I have some information about.

The noise in the data can be modeled as a poisson distribution that blurs signal from a given column in the matrix to adjacent columns. For example, if the original data should be a single peak surrounded by no signal:

0    0    0    1    0    0    0

The measured signal distributed slightly asymmetrically resulting in something like this:

0.001    0.005    0.1    0.5    0.2    0.001    0

If I have a good model of how the noise is distributing the data between the columns, how can I use this information to deconvolve the matrix into an approximation of the original signal?


texts: Filtering voice signal from randomly added frequencies

texts:
I have this assignment that I can't figure how to solve.

I have a sound.wav human-voice file, and use matlab add-random-noise (monotonic)

Goal: to remove the noise from the sound (and the sound from the noise)

What I've done so far is to read the file and then add the noise.
I then do an fft.
Now I need to analyze the frequency spectrum, but I'm not sure how to do it.
And finally do an ifft to recover the sound back.

I've used many kinds of filters: low pass, high pass, band pass, moving average.

The problems were:

I've managed to filter the sound without the noise BUT the sound always got worse as well.
It worked only for specific cases of random noise but not for any random noise that function might produce.

Any help would be appriciated. Thanks!

add_random_noise.m
    function yn=add_random_noise(y,Fs)
        n=[0:length(y)-1]';
        r=rand();
        %%rr=1+rand();
        %%rrr=2+2*rand();
        c=cos(2*pi*440*r*n/Fs);%%+cos(2*pi*440*rr*n/Fs)+cos(2*pi*440*rrr*n/Fs);
        yn=y+c/17;
        yn=yn/max(yn);
    end

include:
add_random_noise.m
sound.wav

EDIT:

I have changed the add "noise" to return only one sinusoid noise.
I run the same noise with my original wav and with zeros and watched the noise in the peak in the range of the lower frequency

now i`m having truble to make the filter:
matlab How to design lpf\bpf\hpf without builtin functions


texts: Why use PSD to analyze a signal instead of its Fourier Transform?

texts:
Let's imagine I have a perfectly known signal $s(t)$ and I want to analyze its frequency components.
I noticed that mostly everyone would use the PSD of $s(t)$ to do that, instead of simply use the magnitude of the Fourier transform of $s(t)$. Is there a reason for that ?


texts: Why remove energy at Nyquist before ifft?

texts: I've heard that you should remove the energy at Nyquist before performing the ifft.  Why is it that you can't leave the Nyquist energy there, is there danger of aliasing somehow?


texts:
I am completely new to image analysis. Do you know how to binarize this image in such a way to get the fibers only?

A step by step procedure would be awesome, but any hint is appreciated as well. I can use Matlab, .NET and I am open to any other tool.

PS: At this address you can find the same question


texts: Real and Imaginary Frequency Responses of a single complex pole

texts:
I'm filtering a real signal with a single complex pole with a complex coefficient (a_re [real part] and a_im [imaginary part]), I also have a gain coefficient but I'm gonna leave it out for the sake of simplicity. The filter has two outputs, one for the real and another for the imaginary part. Here's the filter equation in time for both outputs

Real part output:
$$y_1[n] = x_1[0] + a_{R}  y_1[-1] - a_{I} y_2[-1]$$
Imaginary part output:
$$y_2[n] = a_{I}  y_1[-1] + a_{R} y_2[-1]$$

Where $n$ is the sample number, $x_1[0]$ is the real input, $a_{R}$ & $a_{I}$ are, respectively, the real and imaginary part of the coefficient, $y_1[-1]$ is the previous output of "real part output" and $y_2[-1]$ is the previous output of "imaginary part output".

So, I need to plot the frequency response based on the Z transform. I've already succeeded in dealing with the real version of all this (a real pole), and I'm struggling with the complex version. Anyway, the equation for the real pole is:

$$y[n] = x[0] + a y[-1]$$

Where $a$ is the coefficient of the real pole.

Now, the Z-transform of the real pole is:

$$\tag{1}H(z) = \dfrac{1}{ 1 - a z^{-1}}$$

Where, again, $a$ is the coefficient of the real pole.

In my code I am solving it and dealing  with the z-transform in this way, the input of the function is the $a$ coefficient and the $\omega$ variable (which is the angular frequency in radians per sample). I'm only interested in getting the magnitude output, so I don't care about the phase response of the filter. Anyway, here is how I get it:

$a \rightarrow$ coefficient

$w \rightarrow$ angular frequency

$f_{R} = \cos(\omega)$; $\;\;\;\;f_{I} = \sin(\omega)$

$\mathbf{R} = 1 - a  f_{R}$; $\;\;\mathbf{I} = -a f_{I}$

$\mathbf{Mag} = \sqrt{\mathbf{R}^{2} + \mathbf{I}^{2}}$

$H = \frac{1}{\mathbf{Mag}}$

The whole thing would be:

$$\tag{2}|H(z)| = \dfrac{1}{\sqrt{(1 - a \cos \omega)^2 + (-a \sin \omega)^2)}}$$

Now my trouble is adapting it to the complex version of this z transform, which is supposed to be the same formula as the real pole's z transform in Eq. (1), but I don't know how to adapt the formula/code for the complex pole. Bearing in mind that I want to have the frequency response of the two outputs (the real and imaginary part of the complex filter).

Hopefully, I'd like some help on reaching the frequency response from the coefficient as a function of angular frequency as I've written above for the real pole version of this in Eq.(2)!

Any thoughts, considerations, hints, help is highly welcome.

Thanks a lot


texts: best meathod to lock onto a cdma signal?

texts:
I've been looking at 2 different methods to acquire a known prn code in a CDMA signal. The first trys Xor prn code chip to I sign data. Then average the result over a time period to determine the data bit. I have not had reliable success with this method. The second which I am going to try performs a real multiply of a the code chip and the sign then adds or subtracts the I data magnitude value over a time span.
Is one method known to be more reliable than the other? I'm typing on phone so I can't post any coding right now. If it would help I can update this question later.
Thanks


texts: How is the lowpass to bandpass transformation derived?

texts:
I'm using a lowpass to bandpass filter transformation:

$$ s \leftarrow Q \left( \frac{s}{\omega_0} + \frac{\omega_0}{s} \right) $$

Based on the wiki article: http://en.wikipedia.org/wiki/Prototype_filter which basically is a change of variable. But I was wondering how this transformation is derived?

The reason being I need to create my own transformation so I can combine a low pass with a shelf using a single transformation from a low pass (to use in the Nth order butterworth equation). Knowing how the lowpass to bandpass transformation is derived would take me a step further to experimenting combining other filters.


texts: Normalization purpose in signal processing

texts:
What is the purpose of normalizing the signal?

If we have two signals on hand, how is it used when comparing these two signals?


texts: Is there any way to get intensity value of a certain wavelength for given pixel of a normal  image?

texts:
I understand I can get a color from wavelength but it's not true for the opposite direction.
In case this picture would help.

To make this question simple, let's lighting and calibration aside to an optimized condition.

I wonder if there is any way to get intensity value of a certain spectral wavelength for given pixel of a common image. That is, Wavelength and rgb of a pixel are both known. Based on that, I want to compute intensity of the wavelength for the pixel.


texts: How can I design a PI controller for a closed loop system, which requires a bandwidth equal to the natural frequency and a set Phase Margin?

texts:
For a school project I have to hand in a small report on a this problem. I've been studying so much for it but right now I'm just mixing up everything I've seen and just don't know where to start. I'd really appreciate it if somebody could point me in the right direction.

  For a system with transfer function $$G(s) = 4/(s²+s+4),$$ I have to design a PI controller. The required bandwidth for the closed loop system is the size of the natural frequency. I need to design the PI controller in such way that the phase margin is 50 degrees.


texts: Sliding window minimum and maximum filter in MATLAB

texts: I'm trying to apply a sliding window minimum and maximum filter to an image of a certain window size. Actually, I'm trying to find the optimum window size for it. But I really haven't gotten the hang of it. I presume that I should be using blockproc to implement the sliding window, but not really sure how to find the maximum and minimum filter. As to the implementation itself, should I use loops to slide the window across the entire area of the image ?


texts:
I am trying to understand the OFDM based on DCT instead of DFT in order to compare between their performance. In case of AWGN, the performance is logic, but when I perform channel equalization in DCT-domain, I don't get the performance back ! I think the equalization in DCT domain should be performed in different way. Is that right to perform that process as:

Data sequence -->
S/P -->
IDCT -->
CP insertion -->
P/S -->
Channel + noise -->
S/P -->
CP removal -->
DCT -->
DCT for channel -->
channel equalization -->
received data sequence

That process was performed in Matlab as below :

clear all; clc; clear;

N = 1024;     %number of symbol
M = 2;       %Modulation order
CP = N/4;     %The CP length
Num_sym = 1; %Number of symbols

SNR = 1 : 5 : 35;
BER_OFDM = zeros(1,length(SNR));

for t = 1:length(SNR)
    t
    for t2 = 1:10^1

    L = 32;
    h = randn(L,1) + 1j*randn(L,1);   %Generate the channel

    X = randi([0 1] ,N*Num_sym, 1);        %% Generate data
    X_mod = qammod(X,M);                   %%Data modulation

    X_tr = idct(X_mod);                     %%% IDCT operation
    %%%%
    X_cp = [X_tr(end-CP+1:end); X_tr];       %%Adding CP
    X_chanel = conv(h,X_cp);                 %%Convolution with channel
    X_chanel = awgn(X_chanel,SNR(t),'measured');  %%Add the noise
    Y = X_chanel(CP+1:end-L+1);               %%Remove CP and delay

    H = dct(h,N);                           %Channel in DCT frequency-domain

    EbNo_lin = 10^(SNR(t)/10);                 %linear scale
    Gz = conj(H) ./ (H.*conj(H) + 1./EbNo_lin);   %%Equalizer coefficients

    Y_rec = dct(Y);                            %%DCT operation

    Y_r = Y_rec.*Gz;                          %%Equalization

    sym_r = qamdemod(Y_r,M);                      %Demodularing received signal
    sym_t = qamdemod(X_mod,M);                    %Demodularing transmitting signal

    errors_OFDM = sum(xor(sym_r,sym_t));
    BER_OFDM(t) = errors_OFDM + BER_OFDM(t);

    end

end

err_ber_OFDM = BER_OFDM/(N*log2(M))./10;

semilogy(SNR, err_ber_OFDM, 'k-','LineWidth',1)
hold on
xlabel('SNR (dB)');
ylabel('BER');
grid on

I get such performance as below figure, but I don't think that the performance of OFDM-based DCT is too bad like. I think the channel equalization should be done in different way, but I don't know how


texts: Band Limited Impulse Train - DC offset correction

texts:
I am trying to get band limited saw-tooth wave.
I have found this paper: Alias-Free Digital Synthesis of Classic Analog Waveforms.

I am prototyping this code in common LISP, here is snippet of what I wrote:
https://pastebin.com/Ny4L0EH5

I managed to get something that in time domain looks like sawtooth but it has rising DC offset that causes output to approach infinity (this is generated by blit function).

I read that I should fix this with one pole high-pass filter.
I have put simple IIR but without any good results.
It makes sense for me to remove this DC bias this way since it is digital equivalent of RC differentiating circuit that block DC bias.

I have tried a few different sets of coefficients but without any good results.

I have looked hard with both duckduckgo and google for any information about DC offset correction in this method but i did not managed to find anything.

I have also tried to subtract from output sample average value of one period but again no luck.

I have saved to file band limited version of Dirac comb and tried to filter it with low pass filter in octave but it does not change time domain shape of data in any way, only amplitude decreases. I tried this because I read about this way of performing summation.


texts:
I have a question about changing the phase of a buffer.

Question - How to change the output phase of a buffer by 90 degrees, 180 degrees, and 270 degrees?

I know to change the phase by 180 degrees, we just invert the signal using an inverter. But what about changing it by 90 degrees?

Thank you for your help!


texts: Meaning of Phase response of a filter? In simple words?

texts:
What is meant by phase response of a filter?  Often it is associated with restoration of signal shape,so if a filter has linear phase response, its output looks like input otherwise output is distorted if the filter has non linear phase response

Example of linear phase response filter is FIR filter and IIR filter is example of non linear phase response

But what is meant linear or non linear phase response? Please kindly explain in simple words preferably with an example containing figures or graphs for demonstration


texts: Does changing the frequency of a wave cause a change of the phase of the same wave?

texts:
I know that the phase is the delay between two different waves . But if we change the frequency of two waves by a different value does it cause a varying in the phase between them let say that two waves start from the same point.and if I apply that in the single wave and change the frequency does that related to the phase !
My main question related to the Fm and Pm modulation


texts: What interesting audio can be synthesized with ifft

texts:
I'm working on a presentation on audio synthesis for some fellow software engineers. I want to include some "fourier synthesis" creating time domain samples from frequency domain samples using ifft.

Does anyone know of interesting sounds that can be easily synthesized with ifft?  I can do (bandlimited) square, saw, and triangle waves already, but i'm thinking more along the lines of voice or stringed instruments or something.  Does anyone know of anything that would be a simple to understand but interesting demo for fourier synthesis?

Thanks!


texts: Best Metric to Compare Sparsity of Vectors

texts:
I solved the Basis Pursuit Denoising Problem looking for a sparse solution (I am in compressive sensing):
$$ {x}^{\ast} = \arg \min_{x} \left\{ \frac{1}{2} {\left\| A x - y \right\|}_{2}^{2} + \lambda {\left\| x \right\|}_{1} \right\} $$

for $ 100 $ different $ \lambda $ values with ADMM.

Now I have $100$ different $x^*$ and want to find "the best one" out of them. Moreover, I solved several optimization problems for CS under noise and noisefree cases (e.g. also Lasso, etc.)

What is the best measure to determine sparsity? Found the paper Comparing Measures of Sparsity by Niall Hurley and Scott Rickard, saying it's the Gini index. Do you agree? In my opinion sparsity in this case is not all, isn't it? It's also about solving the equation as good as possible, e.g. having a small norm on:

$$ {\left\| A {x}^{\ast} - y \right\|}_{1} $$

What do you recommend?


texts:
I am trying to solve the following equation for h (an [MxM matrix]):

$$
k[\tau_1,\tau_2]=\sum_{i_1=0}^M \sum_{i_2=0}^M h[i_1,i_2]x[\tau_1-i_1]x[\tau_2-i_2]
$$

I have k, which is a 2D [MxM] symmetric matrix & I have x which is also a 2D [MxM] symmetric toeplitz (autocorrelation) matrix.

I know this is basically a 2D deconvolution problem, but this isnt my field and I cant figure out how to do it in MATLAB.

Also, if possible I would prefer a time domain solution, but frequency domain would also work!

Attempted Solution #1: Division in frequency domain:

h_pred=ifft(fft2(k)./fft2(x));

I think this should work, but for some reason it keeps giving me an imaginary answer!

Attempted Solution #2: Toeplitz matrix inversion:

h_pred=inv(x)*k;

This works better, however while the recovered h_pred is much closer to the true h, it is still very far off...

Please help
Thanks!


texts: Demodulation of a BFSK frequency hopping signal using FFT

texts:
I am not very familiar with DSP techniques, and I am currently working on my research project to develop the firmware of a DSP processor, using C/Linux. It involves the modulation and demodulation of a dataset using Binary FSK and frequency hopping Tx and Rx algorithms. I am able to modulate the signals and have verified the output using MATLAB. However I an a little stuck at demodulation.

The technique that I am using for demodulation is to compare the amplitude of the FFT of the signal received at every hopping sequence and based on that, determine whether the bit received is 1 or 0. I have a predefined function which is used calculate the FFT of the signals. The function is as follows:

void fft(complex_float* in, complex_float* out, int N, int inv)
{
int i, j;
float a, f;
complex_float s, w;
f = inv ? 1.0/N : 1.0;
for (i = 0; i < N; i++) {
s.re = 0;
s.im = 0;
for (j = 0; j < N; j++) {
  a = -2*PI*i*j/N;
  if (inv) a = -a;
  w.re = cos(a);
  w.im = sin(a);
  s.re += in[j].re * w.re - in[j].im * w.im;
  s.im += in[j].im * w.re + in[j].re * w.im;
}
out[i].re = s.re*f;
out[i].im = s.im*f;
}
}

complex_float is a predefined struct :

typedef struct {
float re;
float im;
}complex_float;

In the function fft() N defines the number of FFT points.

The problem that I am facing is how to determine the frequency points of the individual elements in the array returned by fft() function. In MATLAB, this was very easy to do the FFT function in MATLAB took the sampling frequency in account and plotted the power spectrum accordingly. However, In the fft() function defined above there is no frequency involved and hence I am confused as to how to find the power of the signal at any particular frequency, or how to determine the frequency point of any element of the retuned array.


texts: How to find the actual size of an existing object from the pixels of an image?

texts:
I am quite new to vision and OpenCV, so forgive me if this is a stupid question but I have got really confused.

My aim is to detect an object in an image and estimate its actual size. Assume for now I only want length and width, not depth.
Lets say I can detect the object, find its size(length and width) in pixels and I have both the intrinsic and extrinsic parameters of the camera.

The documentation of camera calibration as I understand says that the intrinsic and extrinsic camera parameters can be used to transform from camera coordinate system to world coordinate system. So this should mean converting from pixel coordinates to real coordinates right? And so I should be able to use pixel size and these parameters to find the real size?

But, say the object is photographed at different depth (distance from camera), then its size would come out different using above method. So.....what does it mean transforming from camera to world coordinates?

Any kind of explanation/links/help would be highly appreciated.


texts:
I'm studying about comparison between images to determine quality. I've learned about MSE, SNR and PSNR, and now I'm trying to understand WSNR, which I assume is similar to SNR but with weights.
I have a Python code to calculate the WSNR between two images, but I don't fully understand it and I don't know what is the mathematical definition of this measure. I can't find any paper or book that explain how it is defined and calculated.
Question:
How is WSNR (weighted signal to noise ratio) defined?
Python Code:
import numpy as np
from scipy.ndimage.filters import gaussian_filter as __gaussian_filter
from scipy.ndimage.filters import convolve as __convolve
from scipy.ndimage.filters import correlate as __correlate
from scipy.fftpack import fftshift as __fftshift

def wsnr(reference, query):
"""Computes the Weighted Signal to Noise Ratio (WSNR) metric.

    value = wsnr(reference, query)

    inputs
    ----------
    reference: original image data.
    query    : modified image data to be compared.

    output
    ----------
    value    : wsnr value
    """
    def __genetate_meshgrid(x, y):
        f = lambda u: u / 2 + 0.5 - 1
        (H, W) = map(f, (x, y))
        return (H, W)

    def __create_complex_planes(x, y):
        (H, W) = __genetate_meshgrid(x, y)
        (xplane, yplane) = np.mgrid[-H:H + 1, -W:W + 1]
        return (xplane, yplane)

    def __get_evaluated_contrast_sensivity(plane):
        w = 0.7
        angle = np.angle(plane)
        return ((1.0 - w) / 2.0) * np.cos(4.0 * angle) + (1.0 + w) / 2.0

    def __get_radial_frequency(x, y):
        (xplane, yplane) = __create_complex_planes(x, y)
        nfreq = 60
        plane = (xplane + 1.0j * yplane) / x * 2.0 * nfreq
        s = __get_evaluated_contrast_sensivity(plane)
        radfreq = abs(plane) / s
        return radfreq

    def __generate_CSF(radfreq):
        a = -((0.114 * radfreq) ** 1.1)
        csf = 2.6 * (0.0192 + 0.114 * radfreq) * np.exp(a)
        f = radfreq < 7.8909
        csf[f] = 0.9809
        return csf

    def __weighted_fft_domain(ref, quer, csf):
        err = ref.astype('double') - quer.astype('double')
        err_wt = __fftshift(np.fft.fft2(err)) * csf
        im = np.fft.fft2(ref)
        return (err, err_wt, im)

    def __get_weighted_error_power(err_wt):
        return (err_wt * np.conj(err_wt)).sum()

    def __get_signal_power(im):
        return (im * np.conj(im)).sum()

    def __get_ratio(mss, mse):
        if mse != 0:
            ratio = 10.0 * np.log10(mss / mse)
        else:
            ratio = float("inf")
        return np.real(ratio)

    if not len(reference.shape) < 3:
        reference = __convert_to_luminance(reference)
        query = __convert_to_luminance(query)
    size = reference.shape
    (x, y) = (size[0], size[1])
    radfreq = __get_radial_frequency(x, y)
    csf = __generate_CSF(radfreq)
    (err, err_wt, im) = __weighted_fft_domain(reference, query, csf)
    mse = __get_weighted_error_power(err_wt)
    mss = __get_signal_power(im)
    ratio = __get_ratio(mss, mse)
    return ratio


texts: Machine Learning - any suggestions to solve Python rounding errors?

texts:
I am working with Python to isolate voiced segments from music using the Jamendo corpus for singing detection. Training a model, I break my audio into frames, and have a label (0,1) for each frame.

Unfortunately, due to rounding errors, my labels fall always 1 or 2 frames short.When extracting MFCC features from ~Train Music/03 - sample.ogg, my output is shaped as (13,3709).

This seems reasonable because the length of that file is 3797322 samples and with a hop_length of 1024, the number of frames works out slightly above 3708.

More details about the MFCC feature extraction in the following code snippet:

def MFCC(audio):
    s = []
    for y in audio:
        mfcc = librosa.feature.mfcc(y= y, sr = 16000, n_mfcc=13, n_fft=2048, hop_length = 1024)
        s.append(mfcc)

    return(s)

To now "mark" each frame with the class label (1,0 for singing, no singing respectively), I am converting from seconds to frame numbers but get a maximum value of 3708 (for the longest duration in the file). This is done as per the following function:

def labelToFrames(textCorpus):
    output = []
    for block in textCorpus:
        block_start = int(float(block[0]) * 16000 / 1024)
        block_end = int(float(block[1]) * 16000 / 1024)
        singing = block[2]
        block_range = np.arange(block_start, block_end, 1)  # Step size is 1 (per frame number)
        for x in block_range:
            ms_start = '{0:.3f}'.format(x * 1024 / 16000)           # Converting back to milliseconds
            ms_end = '{0:.3f}'.format((x * 1024 / 16000) + 0.064)   # Converting back to milliseconds
            add_to_output = [ms_start, ms_end, singing]
            output.append(add_to_output)
    return output

I guess librosa's mfcc function is more complex than my 16000/1024 calculation (?). I have tried using Decimal, math.floor and also math.ceil within my block_start and block_stop variables, but I can't seem to match my audio frame length. Are there any other techniques recommended to avoid this error?

EDIT:

To address some discussion points:

y,sr = librosa.load("~Train Music/03 - sample.ogg", 16000, mono = True) followed by y.shape, it yields (3797322,). d = np.abs(librosa.stft(y, n_fft=2048,hop_length = 1024, win_length = 2048))**2 gives a shape of (1025, 3709). Finally performing MFCC on that gives (13, 3709)

Sticking with 44.1khz, the label data is referenced 10 frames short.


texts: Binary robust independent elementary features ("BRIEF") understanding

texts:
How are the descriptors in the BRIEF algorithm matched to each other?

How is the template image found in the another image? How do we compare the descriptors?
I have read that article but didn't understand how they have done it.

They wrote that they compared descriptors by using Hamming distance, but how is invariance to changes in illumination, size, etc. achieved?


texts: How to recover a signal from its complex envelope

texts:
I know there's similar questions here, but I did not fully understand the answers I have read:
in MATLAB I consider the complex envelope of an output signal:

y_low_out   = A_low_out .* exp(1i*Phi_low_in) .* exp(1i*Phi_low_out);

Now, I want to recover a sinusoidal signal, what does y_out look like in terms of superpositions of sine-waves? My guess is:

y_out = A_low_out .* sin(2*pi*f_c*t + Phi_low_in + Phi_low_out)

But this is pure guessing. $f_c$ is the carrier frequency.

Some background to what I am trying to do:
I consider an input signal:

x(t) = A_in .* sin(2*pi*(f_c+f1)*t + Theta1) + A_in .* sin(2*pi*(f_c+f2)*t + Theta2))

This input signal is converted to an equivalent lowpass signal (complex envelope) with amplitude A_low_in and phase Phi_low_in and sent through a nonlinearity which changes both amplitude and phase. The amplitude of my output signal (still complex envelope) is A_low_out, the phase change is Phi_low_out. This description refers to the code that I have written. I have determined A_low_in by taking the absolute value of the analytic signal of x. For some reason A_in and A_low_in are not the same. I figured A_low_in would be twice A_in, but it's not exactly (maybe numerical differences).
Anyways, in the end I get the complex envelope y_low_out. And I want to convert it to y_out.


texts: How do I interpret the DCT step in the MFCC extraction process?

texts:
In most audio processing tasks, one of the most used transformations is MFCC (Mel-frequency cepstral coefficients).

I mostly know the math that's behind the MFCC: I understand both the filterbank step and the Mel frequency scaling.

What I don't get is the DCT (Discrete Cosine Transform) step: What kind of information do I get in this step? What is a good visual representation of this step?


texts: The timing offset remains after Gardner loop

texts:
I would like to use Gardner loop to remove the timing offset in order to demodulate the signal. My signal contains the frequency offset, however I know that the Gardner loop can handle substantial frequency offsets and that is why I apply it.
I have written the following piece of code:
def Gardner_loop(time_support, sig):

    sample_rate = 1/time_support[1]
    sps = int(sample_rate/f_symb) # samples per symbol
    sig_n_1 = int(np.ceil(sps/2)) # x(n-1)
    sig_n_2 = int(sps) # x(n-2)
    loop_filter = f_symb//20 # BW of the loop filter
    loop_num = 10000
    error_symbol = np.zeros(loop_num)

    error_list = []
    sig_wo_time_err = []
    step = 0

    for i in range(1, loop_num):
        error_symbol[i] = sig.real[sig_n_1 + step + int(error_symbol[i-1]*loop_filter)]\
        *(sig.real[sig_n_2 + step + int(error_symbol[i-1]*loop_filter)]\
        - sig.real[step+int(error_symbol[i-1]*loop_filter)])\

        + sig.imag[sig_n_1 + step + int(error_symbol[i-1]*loop_filter)]\
        *(sig.imag[sig_n_2 + step + int(error_symbol[i-1]*loop_filter)]\
        - sig.imag[step+int(error_symbol[i-1]*loop_filter)])

        step+=sps
        error_list.append(error_symbol[i])

        # Update

        sig_wo_time_err.append(sig.real[sig_n_2+step+int(error_symbol[i-1]*loop_filter)] \
        + 1j*sig.imag[sig_n_2+step+int(error_symbol[i-1]*loop_filter)])

    error_to_plot = np.array(error_symbol)
    plt.plot(np.arange(1, loop_num+1), error_to_plot)
    plt.xlabel('loop number')
    plt.ylabel('error')
    return np.array(sig_wo_time_err)

to remove timing offset. The BW of the loop filter I have selected from this article
In this case, sps = Number of samples per symbol = sampling frequency / symbol rate = 4.7. However, it seems to me that the loop is not locked since the error fluctuates as in the picture below:

This is not interpolated signal.
My idea was to increase the number of samples per symbol to reduce the error and lock the loop. After interpolation, the number of samples per symbol = interpolated rate/symbol rate = 5. However, the error still remains and is depicted below:

I have spent many hours trying to tackle on this problem, however I still do not understand what am I doing wrong. So my question is: do I make a mistake somewhere in the code, or do I miss something else to remove the timing offset?
Update:
I have updated the code where I save the signal.
My constellation without interpolation for 1000 samples before timing recovery looks like this:

After symbol timing recovery, again for 1000 points (or for 1000 loops in the Gardner loop)

For the same 1000 loops, my error is represented in the picture below

Note that if I plot more than 1000 samples, constellation becomes not readable anymore (which I think is due to the phase offset)


texts: How can I generate an NxM DCT matrix? In Matlab dctmtx function generates only square matrix (NxN).


texts: Condition for optimal value to be the average

texts: In many single parameter estimation process average is the optimal. Can anyone tell under what condition average is optimal? For example is average optimal if the data is iid, Gaussian or something else?


texts: Erode image until all objects are a single pixel?

texts:
Suppose I would like something more dramatic than skeleton, and would like to erode everything connected to a single point -- but stop there.
As in, I want to keep the single point.  Is this possible with basic morphological operations in OpenCv et al.?


texts: Average Energy of modified QAM Constellation

texts: How to compute the average energy for the modified 16-QAM constellation below? I know that for regular/rectangular 16-QAM we can utilize the 4-PAM along each axes but I'm not sure how can we find $E_{avg}$ here.


texts:
Recently, I have been studying IIR and FIR filters, and trying to create a parametric equalizer using a microcontroller. More specifically, I use the ADC on it to sample the audio at about 44 kHz, and in between grabbing samples, I use the time to process the signal and pass it out on the DAC.
What is a bit confusing to me though is that this process seems to be pretty intense, and often takes most of the time in between samples to do all the processing (with 8 2-pole IIR bandpass filters).
Programs like Spotify, or devices like most Galaxy phones have similar equalizers but they are running on non-dedicated hardware, and I am really confused how they are able to pull this off? Does anybody know if they are basically doing the same thing just faster, or is there some secret sauce to this that I am missing?
I have never tested either of the previously mentioned equalizers in those platforms, so maybe they are just really crappy, but at least just by listening to them, it doesn’t sound horrible.


texts: How to get an objects dimensions from video

texts: I want to be able to measure an object using the camera on my phone. How would I go about getting an objects width and height from a video? My first thought was some type of 3D tracking of the video. But that's complicated and I would like to make sure i'm on the right track. Basically, what would be the easiest way to go about doing this? OpenCV?


texts:
I have solved the problem of finding the filter's order. The order, $N$ is 144 (N= 4/Normalized of BW).

Since the $F_s$(sampling frequency) is 7200Hz, $f_p$(pass-band edge frequency) is 500Hz and $f_s$(stop-band edge frequency) is 700Hz.

Now, I want to:

Create an ideal impulse response i.e, sinc function and then
Apply a Blackman window function to calculate my filter's coefficients.

I am designing the filter in Matlab and I am having difficulty with those two tasks.


texts: Smoothing a discrete function with point wise variances

texts:
In Stark's Introduction to numerical methods [1970] it suggests fitting local polynomials to smooth discrete data before the Fourier transform using least squares minimisation. This gives the following rules for degree three five-point smoothing:

for the first two entries:
$$
g_1 = \frac{1}{70} (69 f_1 + 4 f_2 - 6 f_3 + 4 f_4 - f_5) \\
g_2 = \frac{1}{35} (2 f_1 + 27 f_2 + 12 f_3 - 8 f_4 + 2 f_5) \\
$$

for the last two entries:

$$
g_{N-1} = \frac{1}{35} (2 f_{N-4} - 8 f_{N-3} + 12 f_{N-2} + 27 f_{N-1} + 2 f_{N}) \\
g_{N} = \frac{1}{70} (- f_{N-4} + 4 f_{N-3} - 6 f_{N-2} + 4 f_{N-1} + 69 f_{N})
$$

and for all other entries:

$$
g_{n} = \frac{1}{35} (-3 f_{n-2} + 12 f_{n-1} + 17 f_{n} + 12 f_{n+1} - 3 f_{n+2})
$$

Would a weighted least squares approach following this idea be valid (what would the formulae be?) or are there more widely recognised approaches for smoothing functions with known variances?


texts: How to select the length of a smoothing moving average filter?

texts:
I have a time series signal I need to smoothen to get rid of some noise.

I decided to apply a moving average filter to this signal. I know the choice of the length of the MA is critical.

How should I choose the length?


texts: The normalized cross-correlation of two signals

texts:
I wanted to calculate the normalized cross-correlation function of two signals.
I use the command corr = signal.correlate((s1['Strain']-s1['Strain'].mean())/s1['Strain'].std(), (s2['Strain']-s2['Strain'].mean())/s2['Strain'].std(), mode='full') / len(s1['Strain'])
Both signals have the same length, 30 seconds length and 4096 Hz sampling. s1 and s2 are the GW150914 enter link description here
Here is example data
s1:

            Strain
0        -1.587702e-22
1        -1.425868e-22
2        -1.174897e-22
3        -8.559119e-23
4        -4.949480e-23
.             .
.             .
.             .

I didn't get the normalized cross-correlation, the maximum value is 0.04 in range $[-1;1]$. Where is mistake?


texts:
I want to segment this Red Blood Cell Image on HSI color space.I want to try the Otsu Method but from what I read that method work only on Grayscale image.Do you have any suggestion for suitable image segmentation on those red blood cell image ?
[EDIT]
This is what i want to achieve(Segmented Blood Cell)


texts: How to implement structures for difference equations?

texts:
I want to learn how to design structures that could model difference equations. I don't know how we call this things but i can show a picture:

For example i want to understand the procedure that using addition, multiplication and delay operators like the above i can describe equations like: $ y(k)=1/2(x(k)+x(k-1)) $

Does anyone knows any resources like books, tutorials , articles that explain how to design these structures?


texts: Band-limited random signal with arbitrary distribution?

texts:
I'd like to generate a random discrete-time signal that is band-limited to some bandwidth B (by means of a digital filter, ie in MATLAB). The catch is that I'd like this signal to have an arbitrary PDF, for example, a uniform distribution.

If I start with a uniform distribution, and then filter to band-limit the signal to B, I end up with something that's basically gaussian. Since my band-limiting FIR filter has a large number of taps (say 128), this makes sense given the Central Limit Theorem, as I am basically summing a number of IID random variables.

What I don't understand, is if it is possible to now transform this signal to a new distribution, such as uniform, while maintaining the band-limitedness. I understand that I can transform a RV from one PDF to another by basically integrating/mapping over the CDF (hand-wavy, it's been a while!), but I believe this significantly changes the frequency content of the band-limited signal. For example, if I have a gaussian distribution, and I want to transform to a uniform, I basically map values "near the mean" of the gaussian to something "closer to the max" in the uniform RV. But now the frequency content has changed, and my spectrum will be smeared.

Does what I'm describing make sense, and has anyone looked at something similar? Perhaps what I need to do is look at transforming the signal as a vector, not a scalar (ie treat it as a random process, not individual random variables). But this is where I get a little out of my depth.


texts: Why should the number of subcarriers of OFDM system be a power of 2?

texts:
I notice that almost all papers use the number of subcarriers of 64, 128,256.. so on, which is a number of $2^k$ and $k$ is an integer.

Why is the number of subcarriers always like that? Can we use any number or it must be only using one of those numbers?


texts: Re-recording audio vs. convoluting old audio data - differences?

texts:
I want to transform a studio audio dataset into a smartphone audio dataset (used for a speech recognition application). Would there be any difference in the resulting audio files using one of the following two aproaches?

Replay all audio files and record with smartphone.
Convolute files with smartphone microphone impulse response, assuming the impulse response was determined using the same setup as
    in "1.".


texts:
What is the formal definition of homogeneous noise field? Is there any example of it?

Thanks.


texts: Update Sub Set of the State Vector in Kalman Filter

texts:
I question about Kalman Filter:

If I have system state
$$
\mathbf{X} = [x_1\ x_2\ x_3\ x_4\ x_5]^T,
$$
these state elements are independent. I have measurement from a sensor to correct the predicted state. Normally, the correction will perform on all elements of state $\mathbf{X}$ and produce new covariance matrix $\mathbf{P}$.

If I just want to correct the substate only $x_1$ and $x_2$, and do not do anything with $x_3$, $x_4$ and $x_5$, how can I do that ? and how to compute $\mathbf{P}$?


texts: Chroma-Subsampling: How to properly calculate the data-rate

texts:
I am having a hard time understanding how to calculate the data rate when leveraging chroma-subsampling in e.g. a Y'UV image:

I have the following examples for a calculation:

Image resolution: 352*288
Frequency: 25 fps

For (4:4:4) the example calculation goes as follows:

(352px * 288px) * 3 color channels * 25 fps * 8 bit = 60 825 600 bit/s

So far so good.

But now comes (4:2:0):

(352px*288px) * 1.5 color channels * 25 * 8 = 30 412 800 bit/s

Now, trying to transfer this example to e.g. (4:1:1) I realized I am not sure whether I have a correct understanding how the ratio of 1.5 color channels is being calculated.

My first guess for the calculation was in case of (4:2:0):
2/4*3=1.5 color channels
Likewise for (4:1:1) I would calculate the ratio for the color channels as:

1/4*3=0.75 color channels

But I am simply not sure whether this would be the correct way to go.
As an alternative, I was thinking along following lines:

color channels (4:1:1): 1 Y' + 1/4 UV = 1.25 color channels

Now which would be the correct way of doing this properly?


texts: What is actually transmitted following LPC of a speech frame?

texts:
For each frame, what's sent over to the receiver for decoding? The coefficients , pitch in some bits, voiced/unvoiced classification in another bit?

Another broad question to develop intuition about LPC and CELP for speech encoding.

By using spectral envelope estimation with LPC and related speech coders, does the bitrate reduction happen because the original fluctuations in the frequency spectrum envelope of the original audio signal require more bits to represent? And so in approximating the amplitude spectral envelope fewer bits can be used to represent the speech?

Or, the quality reduction is the result of the elimination of the spectral envelope fluctuations and the bitrate reduction results from only transmitting the LPC coefficients, pitch period, and voiced/unvoiced classification versus transmitting some stream of data for each frame, and then the receiver decodes the coefficients, voiced/unvoiced, and pitch info...?


texts: real refractive index from Kramers Kronig relation

texts:
I have a measurement of the complex part of the refractive index $k$ (where the refractive index is $m = n + i\,k$) measured at a nonlinear grid of wavelengths or frequencies that span several orders of magnitude.

Given $k$, $n$ should be derived from the Kramers-Kronig relations.

However, using a test data set where $n$ is known, I do not recover the correct solution. Any idea what I'm doing wrong (nu = frequency array)?

import scipy.fftpack as ft

plt.loglog(nu, n, label='n')
plt.loglog(nu, k, label='k')
plt.loglog(nu, 1.0 + ft.hilbert(k) / np.pi, label='n from k (orig. grid)')
plt.xlabel('frequency [Hz]')
plt.legend()


texts: I need advice about how to make an audio frequency analyzer

texts:
I just need some clarification of some very basic aspects of creating a frequency analyzer.  I am writing a c++ program, so far it is just a media player, with a custom explorer.  You can make playlists, shuffle, repeat etc.

But my goal is to continue to add functionality.  So far, I just have a signal level meter, and I am now trying to add a frequency analyzer.

I am using fftReal.  So I am at the stage where I have the output, and I am scaling to dbfs like so which I think is correct, but not exactly sure.

for (int i=0; i < fftsize/2; ++i)
    db[i]=20*log(sqrt(imaginary[i]*imaginary[i]+real[i]*real[i]));

So I know I should apply a window function somewhere, and I'm trying to understand how to do that.  But I am wondering this as well.  If I want to display the data as rectangles popping up and down rather than plotting points and connecting with lines, and I want to limit the amount of rectangles to about 50 or so.  Do I take the average db of a range.

I notice that in many applications, the frequency scale is skewed, where the higher frequency, the small distance between units.  What kind of scale is used, and how can I break up the spectrum so that I have a nice visual?

I'm sorry if this has all been discussed before, I have been googling for a while now, and overall I am overwhelmed with information, but there are a just a few key thing which I need to understand first.

Thanks


texts:
I am reading a book on adaptive filters.
For the optimal filter (Wiener), you can analytically obtain the step size for the gradient descent ($0<\mu<\frac{2}{\lambda_{max}}$).
This is not scary, because for the algorithm itself, the correlation matrix is also needed, at the same time it is possible to find the eigenvalues (although in practice it will probably take a long time).
When switching to adaptive filtering (the first considered LMS algorithm), averaging is abandoned (you may not know the statistical characteristics of the process, great!). Now the algorithm uses just instantaneous signal values. But at the same time, the requirements for stability in the book indicate the same. But where can we get them if we have abandoned the knowledge of the correlation matrix, how to correctly implement the algorithm?


texts: Multiply FS Coefficients by 0 or 1 to get Low Pass Filter?

texts:
I have written a bit of code to upsample and interpolate a sample waveform in MATLAB. It's at the point where I have taken a signal, upsampled it 2x and filled in the gaps with 0's, and then I found the fourier series coefficients. I know I need to choose which coefficients to multiply by 0 and which to multiply by 1 in order to get a low-passed signal which I will then perform an inverse Fourier Transform on, and this will give me back my interpolated upsampled signal.

First of all - why do I need to low pass this signal? I suppose I kind of added period information my bringing in those 0's and upsampling the signal, but I don't understand how to find a cutoff for this ideal lowpass I need, and how to translate that cutoff frequency into 0's and 1's scaling my FS coefficients. I heard I need to multiply every other sample by 0, but I don't know why this would be correct.

Here is my code:

function [a] = UpSampler(x)
N=length(x);
Nu = 2*N;
xt=zeros(1,Nu);

for n=1:Nu
    if(mod(n-1,2)==0)
        xt(n)=x(1+(n-1)/2);
    end
    if(mod(n-1,2)==1)
        xt(n)=0;
    end
end

a = getDTFS(xt);

%multiply a coefficients by 0 or 1
%Perform Inverse FT to get interpolated signal

end


texts: How to derive the bin width / maximum frequency of DFT/STFT

texts:
The bin width is given by the $F_s / N$ and the maximum frequency is given by $F_s / 2$ where $F_s$ is the sample rate and $N$ is the number of samples of the DFT.

Going between these is trivial, but I'm having trouble seeing why these formulas are correct.

I see that these formulas imply that for this image:

  halving the window duration would remove every second frequency bin, but keeping the same window duration and removing every second sample (halving the sample rate) would correspond to removing the right half of the bins.

Edit

Thank's to Dan's answer I now understand that the DTF is $N$ bins periodic ($N/2$ bins periodic for a real signal). But I want to see why either the bin width is $1/T$ where $T$ is the duration the samples are taken over, or why the DTF is $F_s$ Hz periodic. (I think bin width as $F_s / N$ is a red herring)

Aside
Sam informed me that my question assumes that the DFT size is equal to the window length and this is just a special case.


texts: How to deal with the mirrored components of a FFT? And another question

texts:
I'm trying to simulate the behavior of my FFT spectrum analyzer, I am using Mathematica, but my questions are more conceptual than technical implementation of code...I hope.

If we consider a simple cosine wave transient signal $y(t) = A \cos(2 \pi \nu_0 t)$, where $A = 1$ and $\nu_{0} = 24$ $\rm{Hz}$. This signal is sampled at some time interval $dt$ and this gives me 2048 points (2048 because of my device), I want to resolve this data to some frequency resolution, $\Delta f$, and from the Fourier limit: $t_{sample} = 1/\Delta{f}$ so my sample rate, $S_{r} = 1/dt$ is adjusted accordingly.

I then take the absolute value of the FFT of these 2048 points and reconstruct the frequency component using $\Delta f$ (in this example $\Delta f = 0.0625$ $\rm{Hz}$, this value is because of my device and the features I want to resolve) and the index of the resultant FFT, and get this:

I understand that the symmetric appearance of an FFT comes from the real and imaginary components having the same response to signals -- the absolute value of their responses will be the same. But the phase of the real and imaginary parts is shifted by $\pi / 4$. This gives us our mirrored appearance.

FIRST QUESTION: Usually I would just throw the right hand frequencies away and multiply the amplitude of what is left by $2$. Is this the correct approach or should I reverse the order of one half of these 2048 points and then add them together (intuitively I would say no, because both halves should be identical so multiplying by $2$ should suffice)? What is the best method of reconstructing the frequency component/axis?

SECOND QUESTION: This is something I just don't understand. If I increase the frequency of my transient signal $\nu_{0}$ then the two mirrored peaks move closer to each other, converging in the center. This makes sense. If I keep increasing $\nu_{0}$ of the transient arbitrarily (past $124$ $\rm{Hz}$ which corresponds to the limit of frequencies) however I still see a peak, in fact I see this oscillating behavior of the peaks just bouncing back and forth as I further increase the frequency -- Why?! Is this natural or because of the way I define my frequencies.


texts: Wavelet packet analysis, Sampling rate, Decomposition level

texts:
I am running Wavelet packet analysis for a 1-D signal. I want to estimate dominate frequencies in the original signal.
Imagine I have following signal:

fn=[0.4 0.6 1.0];E=[0.03 0.05 0.1];A=[1.0 2.0 3.0];
wn=2*pi.*fn;a=E.*wn;wd=wn.*sqrt(1-E.^2);
y=A(1)*exp(-a(1)*t).*sin(wd(1)*t)+A(2)*exp(-a(2)*t).*...
 sin(wd(2)*t)+A(3)*exp(-a(3)*t).*sin(wd(3)*t);

so it contains these frequencies  $\big\{0.4Hz,0.6Hz, 1Hz\big\}$
My question is what is the best level of decomposition to see these frequencies?
I select $10Hz$ for sampling rate.
So my question, what sampling rate should I select for this case to identify these 3 freq?
(I've select 10Hz for sampling rate)
and what level of decomposition should I select to see all of these 3 frequencies?
Here is my little scrip for this signal with wavelet packet

df=30;N=300;                            //df:final time,N:segments
t=linspace(0,df,N); dt=t(3)-t(2);      // dt: time intervals
fs=ceil(inv(t(2)-t(1)));              // fs: sampling freq
fn=[0.4 0.6 1.0];E=[0.03 0.05 0.1];   //fn:freq, E:damping
 A=[1.0 2.0 3.0];wn=2*pi.*fn;a=E.*wn;wd=wn.*sqrt(1-E.^2);  // A:Amplitude.
y=A(1)*exp(-a(1)*t).*sin(wd(1)*t)+A(2)*exp(-a(2)*t).*...
      sin(wd(2)*t)+A(3)*exp(-a(3)*t).*sin(wd(3)*t);
level = 3;
wpt = wpdec(y,level,'coif5');
[Spec,Time,Freq] = wpspectrum(wpt,fs,'plot');

here are result for Level 3,4,5,6 decompositions.

 as you see even in level 6 it does not show accurate frequencies. ans If I increase the level to 10 I will observe higher freq that are not exist in the actual signal.
for level 10:


texts: How to reduce reconstruction noise from a filtered back-projection reconstruction of a circularly symmetric image?

texts:
I need to perform tomographic reconstruction of an axially symmetric object. Because of the axial symmetry only one projection of the object (one angle) was taken. I implemented the filtered back-projection algorithm in MATLAB. However, when I run this algorithm on a synthetic image which has circular symmetry, but with added noise to its projection, I get a lot of noise in the reconstruction, especially in the center of the reconstruction.

Is there a way to modify the filtered back-projection algorithm so that the noise in the reconstruction will be reduced while this reduction will not result in loss of information ??  or is the filtered back-projection algorithm just doesn't work for a noisy image which has circular symmetry and as a result the reconstruction is performed with only one projection ??  (I must add that I added the noise to the synthetic image projection in a symmetric image).


texts: Frequency-domain zero padding - special treatment of X[N/2]

texts:
Suppose we wish to interpolate a periodic signal with an even number of samples (e.g. N=8) by zero-padding in the frequency domain.

Let the DFT X=[A,B,C,D,E,F,G,H]
Now let's pad it to 16 samples to give Y.  Every textbook example and online tutorial I have seen inserts zeros at [Y4...Y11] giving
Y=[2A,2B,2C,2D,0,0,0,0,0,0,0,0,2E,2F,2G,2H].
(Then y = idft(Y) is the interpolated signal.)

Why not instead use
Y=[2A,2B,2C,2D,E,0,0,0,0,0,0,0,E,2F,2G,2H]?

As far as I can tell (my math knowledge is limited):

It minimises the total power
It ensures that if x is real-valued then so is y
y still intersects x at all the sample points, as required (I think this is true for any p where Y=[2A,2B,2C,2D,pE,0,0,0,0,0,0,0,(2-p)E,2F,2G,2H])

So why is it never done this way?

Edit: x is not necessarily real-valued or band-limited.


texts: Power spectral density: where does the difference from decibel and regular representation come from?

texts:
I am generating a simple DC signal with noise. As expected the fft shows that there is one peak at 0Hz and some noise. I am also trying to get the spectral power density representation of that signal, and I came to a problem.

When i am using a logarithmic scale, to show power to frequency relation in dB/Hz, I am observing a gradually decreasing curve (second curve in the picture below). It shows that the lower the frequencies the less power they have in relation to same value. I don't understand it  - there is supposed to be only one peak at 0Hz, all other frequencies should have the same power.

This phenomenon can be observed on a regular graph, on non-logarithmic scale. This other diagram (first curve) shows clearly that there is just one frequency with "energy".
Where does this difference come from?
Does it come from leakage? If so, then why it is not relevant for the non-logarithmic spectrum?

The matlab code below generates the signal and the graphs.

clear
home
close all

signal_multiplier = 3;
noise_bits =8;
signal_bits = 11;

fs = 1;
T = 1/fs;
t = 1:T:(2^signal_bits)-T;
noise_vec2 = randi([-(2^noise_bits) 2^noise_bits],1,length(t));

signal = signal_multiplier*t + noise_vec2;

NFFT = length(signal);
[P, F] = periodogram(signal,[],NFFT,fs,'power');
PdBW = 10*log10(P);

figure(1)
subplot(2,1,1);
plot(F,P)
grid on
title('PSD')
xlabel('Frequency (Hz)')
ylabel('Power/Frequency')
subplot(2,1,2)
plot(F,PdBW)
grid on
title('PSD using logarithmic values()')
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB/Hz)')


texts: What is the difference between FIR1 and FILTFILt in MATLAB?

texts: In MATLAB there are filters with 0 phase filtfilt and ordinary FIR1. When I filtered a group of simulated data, I found that filtFilt was superior to FIR1 in terms of phase delay. I know the difference. But would there be any essential difference between filtfilt-designed filter and FIR1-designed filter when written into FPGA? Or are they just different in Matlab when it is used to filter data?  Besides, could you please tell me how to set the order of filtfilt? I would like to compare the filtering performance difference between FILtFilt and FIR1 of the same order. Thanks!


texts:
I want to implement the following algorithm. Let  $I$  denote  noisy  image  and  for  each  pixel of  $I$  is represented by $I(x,y)$. A Sliding window of size $3\times3$ centered at $I(x, y)$ is defined. At first we estimate that the center pixel of window is noise or not. To differ the noise and signal, we calculate the mean ($m$) and standard deviation ($s$) of the window. If the pixel is between $m-s$ and $m+s$ then it is a signal otherwise noise and we will apply weiner Filter.

Algorithm:

1.  Read an Image I.
2.  Take a Sliding Window or Mask of size 3X3.
3.  Calculate the Mean (m) and Standard deviation (s) from the Mask.
4.  Calculate the threshold as follows
         t1 = m-s and t2 = m+s

5.  IF     t1 <= I(x,y) and  I(x,y)  <= t2 THEN

           Result(x,y) = I(x,y)
    ELSE
        Result(x, y) = weiner2(I(x,y)

6.  Repeat the step 3 , 4 and 5 on Entire Image.
7.  Display the Resultant Image.

First, I tried to implement this algorithm for median filter, so that I can extend it to other filters such as wiener, anisotropic diffusion,butterworth etc.  I did implement it correctly with the help of 'tbirdal'. Thank you 'tbirdal'. But I am unable to extend the implementation for other filters.

**I don't want to use loop.. Because it takes more time to run .. 20 mins to run a $512 \times 512$ image.. I want to implement this using functions such as nlfilter, imfilter etc **
Program:

A=imread('npan3.jpg');

fun = @(x) noisedetect_2(x(:));
B = nlfilter(A,[3 3],fun);

imshow(A), figure, imshow(B, [])

The noisedetect_2 function is given as

function nd=noisedetect_2(sub)
sub=vec2mat(sub,3);
[row col]=size(sub);
%siz=size(sub)
s=std2(sub);
m=mean2(sub);
t1=m-s;
t2=m+s;
nd=zeros(row,col);

if (sub(2,2)>=t1 && sub(2,2)<=t2)
    nd=sub;
else
    nd=wiener2(sub);
end

I get the following error:

Assignment has more non-singleton rhs dimensions than non-singleton subscripts

Error in nlfilter (line 75)
        b(i,j) = feval(fun,x,params{:});

Error in nd_test_2 (line 10)
B = nlfilter(A,[3 3],fun);

I understand what the error means. the nlfilter is expecting single value but the function returns $3\times3$ matrix. How to solve this issue?. I just want to apply any filter only when the pixel falls out of the range $t_1$ and $t_2$ because it is considered as noisy pixel. Otherwise, it is retained as such because it is considered as noisefree pixel.

Please help me implement the algorithm.


texts: Finding the coefficients of the block diagram

texts:
I was looking at the question from a textbook of mine and it doesn't make too much sense to me:

Design a ﬁrst-order ﬁlter, linearly extrapolating the current and
previous measurement. This is equivalent to the assumption that the
upcoming measurement will be on a straight line connecting the current
and previous measurement. The output of the ﬁlter shall anticipate
this by giving this output now.

The coefficients (2 and -1) are given in the answer, but I still don't understand how they were found.
I know that the general equation of a digital filter is given as:
\begin{equation}
    y[k] = \sum_{i=0}^{N} b_i u[k-i] - \sum_{j=1}^{N} a_i y[k-j]
\end{equation}
How does one go from that, using the description of the problem, to:
\begin{equation}
    y[k] = 2 x[k] - x[k-1]
\end{equation}
This is the answer given in the book:


texts: Separating 'signal' from periodic 'noise'

texts:
I have data in which the signal we want is a few gaussian peaks superposed on a noisy background in which some of the 'noise' is periodic.  Simply suppressing the associated fourier modes is effective (e.g., a 'notch filter') but also suppresses the signal, since it has some power associated with those frequencies.

What is a good approach for removing these periodic signals without removing the corresponding signal?

This figure shows the data (top; with X-axis labeled by pixel #) and the absolute value of its fourier transform (bottom; X-axis is pixel frequency) with relevant regions identified.

An approach I've considered is least-squares fitting a series of sine functions with frequencies limited to the identified range and the signal-containing region masked out, but this is computationally expensive and it is difficult to automatically identify the affected frequency range.


texts: Causality, Region of Convergence and which side is the system on

texts:
I want to try to summarize some relationships between z transform, ROC and causality. But I am getting confused between left/right sided sequences, the ROC and causal.

If you need stability then the ROC must contain the unit circle. I can understand that.
If you need a causal system then the ROC must contain infinity and the system function will be a right-sided sequence. I am fine with  that too.
If you need an anticausal system then the ROC must contain the origin and the system function will be a left-sided sequence. Yes that's fine.
If you need both, stability and causality, all the poles of the system function must be inside the unit circle.

But #4, doesn't it violate number 2 because we don't include infinity?

Got it from here: http://en.wikipedia.org/wiki/Z-transform#Example_2_.28causal_ROC.29

Thanks


texts:
This article here talks about a room impulse response generator which takes the source and destination coordinates, sound speed, sampling rate, room dimension and wall reflection coefficients as input and outputs the impulse response from source to destination.

The equation (36) for the same calculation is given below.

My question is if we know different room impulse responses for different source-destination combinations and other constant parameters like room dimension, sound speed etc. except the wall reflection coefficients($\beta_1,\beta_2,..,\beta_6$). Then is it possible to find the $\beta$ values from the Room Impulse Responses(RIRs) and source destination coordinates?.

What I am essentially trying to do is measure some real RIR of a given room of which dimension and source-destination coordinates are known to me. I want to build a model out this real room measurements with minimal model  mismatch. How can I do the same?.

Thanks & Regards,


texts: Unable to estimate for AR model using OLS, Yule Walker and MLE

texts:
I am  learning estimation methods following the book of Steven Kay, "Fundamentals of Statistical Processing, Volume I: Estimation Theory "
Theory says that if the measurement noise is Gaussian, the OLS = MLE. However, in the noise-free case, why am I not getting the true estimates? I have not checked for the noisy case. Can somebody please help in correcting my code and the concept as to why none of the estimation methods such as the yule-walker, OLS and MLE are working to give the true estimates of an AR(2) model. Where am I going wrong?
clear
h =[1,0.195,- 0.95]; %true coeffficients
L=length(h);
N=256;
 x=round(randn(1,N)); %driving signal

y  = filter(1,h,x); % generate the AR(2) model
arcoeffs_yule = aryule(y ,2)
arcoeff_mle = mle(y)
Y=y';
x=x';
arcoeffs_ols =  [Y -[0;Y(1:end-1)] -[0;x(1:end-1)]]\x

arcoeffs_yule =

    1.0000    0.9285    0.0000

arcoeff_mle =

   1.0e+07 *

   -0.0941    7.8185

arcoeffs_ols =

    0.1405
   -0.1513
    0.1527

UPDATE to the code: based on the comment, I have put the length of the AR model=3 inside aryule instead of passing the order 2. Then I have passed some initial values to the filtered signal. Yet getting incorrect estimates.
   clear
    h =[1,0.195,- 0.95]; %true coeffficients
    L=length(h);
    N=256;
     x=round(randn(1,N)); %driving signal

y(1) = 0.1;
y(2) = 0.2;

% Generate the AR model.
for i =3 : N
    y(i) = 0.195 *y(i-1) -0.95*y(i-2) + x(i);
end

%y = filter(1,h,x);

arcoeffs_yule = aryule(y,L)
arcoeffs_mle = mle(y)
Y=y';
x=x';
arcoeffs_ols =  [Y -[0;Y(1:end-1)] -[0;x(1:end-1)]]\x

arcoeffs_yule =

    1.0000   -0.0420    0.9202    0.1261

arcoeffs_mle =

    0.0040    2.9161

arcoeffs_ols =
  0.1108    0.0263    0.1613


texts:
What I am attempting to do is take a signal (.wav audio file) which contains calls from different Bat species. Here is a graph representation of the signal:

What I have done so far is:

1) Emphasise the signal

2) Compress the signal

This is the output I got plotted:

I need to remove the signals that do not warrent consideration. At the moment, the signal is in a 1D representation, and, after the Separation has been completed the signal will be a 2D representation containing JUST the calls, thus I can perform analysis on each of the calls.

This is my suggested methodology to do this:

1) Split the signal into blocks of: 256x100 (so they overlap)

2) Calculate the RMS (Root-mean-square) of each of the blocks

3) Compare (2) to a threshold value, and, if the block warrants consideration then this can be classed as a call.

But the problem is this:

1) Will splitting the signal into blocks of 256x100 go over the calls and thus making some of the calls lost when calculating the RMS?

2) How do I therefore pick a suitable threshold value? Would finding the peaks in the signal be an accurate assumption for this?

I hope someone can help me :)


texts: Re-implementing the "spectrogram" function from matlab

texts:
I am trying to make a spectrogram viewer without using the spectrogram command. For this purpose, I used a sin function. I broke up the input signal into 256 segments, multiplied each segment with the Hamming window, took the FFT of each slice and then stored them to obtain spectra of the signal.

My problem is that I couldn't get the exact spectrum. It is probably due to the fact that I could not divide input signal into overlapping segments?

This is my code: Can someone help me where I am doing it wrong?

Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = 10240;                     % Length of signal
NFFT=10240;                   % number of fft points
t = (0:L-1)*T;                % Time vector
x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
X = hann(L).*x';
figure
plot(X)
k=256;

for idx = 1:length(x)/k
     %take a slice from the input vector
    slice =  x(1+(idx-1)*k:idx*k);
    %multiply it with the window and transform it into frequency domain
    spectrum =  fft(slice'.*hamming(k),k);
    %get the spectrum magnitude at each of the 256 frequency points and store it
    mag_spectrum(:,idx*1024) = abs(spectrum).^2;

end

figure
w=hamming(k);
plot(w)
figure
plot(t,X,'k');
xlabel('Time (sec)','FontSize',14);
ylabel('\it{x(t)}','FontSize',14);

% spectrogram(x,128);
% SET PARAMETERS

figure
plot(slice); % To look for it is windowed or not which is not in my case?

%figure
%plot(mag_spectrum);
figure
spectrogram(X,256); % Real spectogram


texts: Difference equation with variable coefficients in MATLAB

texts:
Given the following difference equation

$$y[n]-\left(\dfrac{n}{n+1}\right)y[n-1] = n\cdot x[n]$$

How can we use MATLAB to solve it?

I know if the coefficients are constant we can simply use filter(b, a,x), but not sure how to do it with variable coefficients.


texts: Dealing with bias in an FFT for frequency correction

texts:
I'm working on using an FFT for frequency correction of a generic modulation signal. My algorithm basically attempts to balance the energy in the positive and negative frequencies of the FFT to be equal. One thing I'm wondering is if I have an inherent bias. Noting the following ordering returned by the bins of the FFT with N even:

After computing the FFT, my algorithm splits up the positive side into the bins $X_{up} = X[0] +...+ X[N/2 - 1]$ and the negative side into the bins $X_{low} = X[N/2] +...+ X[N-1]$ and then form an error signal as $e = X_{up} - X_{low}$
However, I realized this might be resulting in some bias because my positive frequencies are actually including DC, while my negative frequencies are including the nyquist component. So I'm wondering if the correct solution here is to drop both the DC and nyquist bins (I could add them to both the upper and lower terms but since I substract them they'd just cancel anyways) and instead do the following:
$X_{up} = X[1] + ... + X[N/2 - 1]$ and $X_{low} = X[N/2 + 1] +...+ x[N-1]$. This results in an equal number of bins, N/2-1, for both the lower and upper terms and I assume will be unbiased?


texts: Why is histogram equalization called as "equalization"?

texts:
When output histogram is not same as reference histogram, why is the method called histogram equalization?

For example below are two figures from http://en.wikipedia.org/wiki/Histogram_equalization. The output histogram (2nd figure) looks nothing like a discrete uniform distribution. Why is it called "equalization" then?