texts: Edge artifacts in frequency-domain Moving Average filter

texts:
I am currently trying to implement in python an algorithm to identify and filter out one specific periodic EEG artifact -- the grad fMRI artifact -- as described in this paper, and I am incurring in what I believe to be strong edge artifacts with little to no clue on how to take care of them.
As far as my understanding, the basic idea described there is an improved version of the Average Artifact Subtraction (AAS) based on the idea of Optimised Moving-Average (OMA) that can be briefly outlined as follows.
The process finding the best artifact template is performed through iterations, and the Fourier transform of one single step of the iteration is:
$$
1 - H(z) = 1 - \frac{1}{M^2} \cdot \frac{(1-z^{-M})\cdot(1-z^M)}{(1-z^{-1})\cdot(1-z)}
$$
M being the period of the artifact.
The transfer function for the filter implying J iterations of the moving average then becomes:
$$
H_J = 1- (1-H(z))^J
$$
Finally, in order to obtain a more refined filter, the above filter is applied in a series of J cascades leading to the final transfer function:
$$
H_L(z) = (H_J(z))^L
$$

The computational implementation of the filter has been approached in the frequency domain and tested on a data from a measurement of about $1.8 \cdot 10^2$ s, where my sampling frequency was of $5 \cdot 10^3$ Hz and the artifact had a period of roughly $2\cdot10^{-2}$ s.
To avoid violations of the stationarity of the signal, the original signal was high passed to make sure it has zero mean.
The time domain is then transformed into the frequency domain as follows:
N = len(signal)
k = np.arange(N)
z = np.exp(1j * 2 * np.pi * k / N)

and the signal transformed
signal_tilde = fft(signal)
and the result is taken to be the anti-transform of the multiplication of the transfer function with the signal (in frequency space).
signal_filtered = ifft(H_J(z)*signal_tilde)
Without any padding, the resulting signals ends up with a lot of edge artifact.
In order to try to avoid it, the signal was padded  with symmetric padding before starting any
np.pad(signal, (pad_before, pad_after), mode='symmetric')
And consequently the outcome signal has been taken to be the cropped
signal_filtered = np.real(signal_filtered[pad_before:pad_after])

But again the transformed+cropped signal shows unequivocal signs of strong artifacts on the edges covering approximately 3% of my whole signal. One single channel is shown in the figures below for clarity.

Does anyone have any further suggestion or hints on what I could try to get rid of edge artifacts in this case?


texts: How do I retrieve texture using GLCM and classify using SVM Classifier?

texts:
I'm on a project of liver tumor segmentation and classification. I used Region Growing and FCM for liver and tumor segmentation respectively. Then, I used Gray Level Co-occurence matrix for texture feature extraction. I have to use Support Vector Machine for Classification. But I don't know how to normalize the feature vectors so that I can give it as an input to the SVM. Can anyone tell how to program it in Matlab?

To the GLCM program, I gave the tumor segmented image as input. Was I correct? If so, I think, then, my output will also be correct.

My glcm coding, as far as I have tried is,

I = imread('fzliver3.jpg');
GLCM = graycomatrix(I,'Offset',[2 0;0 2]);
stats = graycoprops(GLCM,'all')
t1= struct2array(stats)

I2 = imread('fzliver4.jpg');
GLCM2 = graycomatrix(I2,'Offset',[2 0;0 2]);
stats2 = graycoprops(GLCM2,'all')
t2= struct2array(stats2)

I3 = imread('fzliver5.jpg');
GLCM3 = graycomatrix(I3,'Offset',[2 0;0 2]);
stats3 = graycoprops(GLCM3,'all')
t3= struct2array(stats3)

t=[t1,t2,t3]
xmin = min(t); xmax = max(t);
scale = xmax-xmin;
tf=(x-xmin)/scale

Was this a correct implementation? Also, I get an error at the last line.

My output is:

stats =

   Contrast: [0.0510 0.0503]
Correlation: [0.9513 0.9519]
     Energy: [0.8988 0.8988]
Homogeneity: [0.9930 0.9935]
t1 =

Columns 1 through 6

0.0510    0.0503    0.9513    0.9519    0.8988    0.8988
Columns 7 through 8

0.9930    0.9935
stats2 =

   Contrast: [0.0345 0.0339]
Correlation: [0.8223 0.8255]
     Energy: [0.9616 0.9617]
Homogeneity: [0.9957 0.9957]
t2 =

Columns 1 through 6

0.0345    0.0339    0.8223    0.8255    0.9616    0.9617
Columns 7 through 8

0.9957    0.9957
stats3 =

   Contrast: [0.0230 0.0246]
Correlation: [0.7450 0.7270]
     Energy: [0.9815 0.9813]
Homogeneity: [0.9971 0.9970]
t3 =

Columns 1 through 6

0.0230    0.0246    0.7450    0.7270    0.9815    0.9813
Columns 7 through 8

0.9971    0.9970

t =

Columns 1 through 6

0.0510    0.0503    0.9513    0.9519    0.8988    0.8988

Columns 7 through 12

0.9930    0.9935    0.0345    0.0339    0.8223    0.8255

Columns 13 through 18

0.9616    0.9617    0.9957    0.9957    0.0230    0.0246

Columns 19 through 24

0.7450    0.7270    0.9815    0.9813    0.9971    0.9970

??? Error using ==> minus
    Matrix dimensions must agree.

The input images are:


texts:
I have recently started to use MathNet.Numerics, and one of my needs is to perform FFT analysis of Electromyographical signals.

The first test I did was to get some 15 seconds of signals at 2kHz, this is how it plots:

I have a Python script I wrote some time ago, and the FFT spectrum generated is below:

lenfft = len(signal)/2.
nyq = 2000/2.
freq_axis = numpy.linspace(0, nyq, lenfft+1)
spectrum = numpy.abs(numpy.fft.rfft(signal))
plt.plot(freq_axis, spectrum)
plt.show()

The problem is, when I use MathNet.Numerics, I get a weird result (code and image below):

Complex[] signal = ArrayOfDoublesFromFile()
                        .Select(val => new Complex(val, 0))
                        .ToArray();

Fourier.Forward(signal, FourierOptions.Default); // inplace

SaveToFile(signal);

I know that if the inputs are pure real numbers, the output should be symmetrical or something like that, but since Fourier.Forward() requires a MatNet.Numerics.Complex[] array, I am creating an array with imaginary part set to zero, and also using only real part of result for plotting, is that correct?

So, the question is: whats the meaning of the output generated by MathNet.Numerics.Fourier.Forward(), and what should I make if I want a result in the visual form most commonly associated with FFT, that is, the middle figure?


texts: Converting complex baseband signal to passband signal

texts:
I have to convert baseband signal (BW-20MHz) to passband signal centred at 20MHz.
I am implementing this in matlab. so, this is how I tried it.
I have already generated LTE Waveform BW- 20MHz , fs =30.72e6, size(LTE Waveform)= 307200 samples.

n = [0:size(LTE_Waveform)-1]';
x = exp(1i*2*pi*20e6*(1/30.72e6)*n);

z =LTE_Waveform.*x;

fvtool(z, 'fs',30.72e6)

This is the LTE waveform

and the passband spectrum I got looks like this

How can I get the spectrum properly centred around 20MHz?


texts: How to identify stable deltas in a noisy signal

texts:
This is a problem I have had to deal with in several different contexts now.  Given a real time signal which you expect to act as a noisy step function (periods of stability separated by quick changes in magnitude).  How does one identify when the signal has stabilized?  The tension in this problem is between identifying stability as fast as possible and reporting a new stable value with as much accuracy as possible.

My solutions have varied depending on which of those I felt was more important and on how good my transition model was.  I have taken to calling this the stable deltas problem but have been unable to find any references on how to approach it.  Is this a studied problem in signal processing?  Does it have a name?


texts: Translation invariance in max-pooling and cascading with convolutional layer

texts:
I am reading through this page about convolutional neural networks and I am confused about the part about max-pooling (downsampling a signal/image with the maximum value in a block).

Apparently, one of the advantages of max-pooling is translation invariance:

  It provides a form of translation invariance. Imagine cascading a max-pooling layer with a convolutional layer. There are 8 directions in which one can translate the input image by a single pixel. If max-pooling is done over a 2x2 region, 3 out of these 8 possible configurations will produce exactly the same output at the convolutional layer. For max-pooling over a 3x3 window, this jumps to 5/8.

I kept thinking about this and I can't figure out why this is true.

I am assuming the author means doing a translation to the left/right of the image, passing through a max-pooling layer and convolving. Then comparing the outputs of the convolution with the same process but for a different translation direction

I even did a Python script to illustrate the process of max-pooling -> convolution with 8 different translations by one pixel and the original matrix

from itertools import product

from cv2 import warpAffine
from matplotlib.pyplot import figure, gca, imshow, show, subplot
from numpy import asarray, concatenate, eye, float32, max, min, mod, reshape, zeros
from numpy.random import binomial
from scipy.signal import convolve2d
from skimage.measure import block_reduce

def list_to_2d(x):
    return reshape(asarray(x), (len(x), 1))

def npix_to_m(x):
    return concatenate((eye(2), list_to_2d(item)), axis=1)

def show_images(xl):
    vmin = min(asarray(xl))
    vmax = max(asarray(xl))
    figure(figsize=(4.6,4.6))
    for (x, i) in zip(xl, range(len(xl))):
        subplot(3, 3, i+1)
        gca().xaxis.set_major_locator(plt.NullLocator())
        gca().yaxis.set_major_locator(plt.NullLocator())
        imshow(x, interpolation='None', cmap='gray', vmin=vmin, vmax=vmax)
    show()

x = float32(binomial(3, 0.5, (10,10)))

translations = list(product((-1, 0, 1), repeat=2))

M_list = [npix_to_m(item) for item in translations]

x_trans = [warpAffine(x, M, (10,10)) for M in M_list]
show_images(x_trans)
print 'Original matrices\n'

x_maxp = [block_reduce(x, block_size=(2,2), func=max) for x in x_trans]
show_images(x_maxp)
print 'Max-pooled matrices\n'

filt = float32(binomial(3, 0.5, (2,2)))
x_conv = [convolve2d(x, filt, mode='valid') for x in x_maxp]
show_images(x_conv)
print 'Convolved max-pooled matrices\n'

These are the outputs:

As you can see none of them are exactly the same. So what am I not understanding here?


texts:
I have this problem where I have multiple points in a 2D space. A subset of these points makes up a rectangular grid. Some points of the grid might be missing. The shape and number of points the grid consists of are known. Also, the other points not making up the grid are few. The general location and rotation of the grid are also known. But we don't know the exact location or rotation. What would be a good, and most important, efficient way to find this grid?

A perfect grid:

The grid can appear as shown below with a missing point:

Or with for example a rotation and being split in two as shown below:

In the actual problem, the grid is larger but always rectangular, also we can assume the gird is located around the center of the data.

The points come from an image where there is monospaced text. The missing points and noise are due to glare in this image. The points are the centers of the
(by preprocessing) found characters (so not exact). The text in the image always consists of the same number of characters and lines. The idea is reconstructing this grid to be able to find the lost characters. We do not know what points are characters and what points are not, but I am aware this could be approximated by analyzing the "blobs" we find in the image.


texts: What is the intuition behind convolution?

texts: I have been using convolution for finding outputs of various systems .I know how to use it.But I still don't know what does convolution exactly  means? how one can define convolution ?


texts:
The title might be unclear, but the problem is this. I have a signal sampled 1500 times with a rate of 60/s, and a sensor array 512 units large. There is a lot of noise, echo and other frequencies being picked up, but I am interested in only one. First I do a spike removal, then a bandpassfilter (butterworth) around the frequency range I suspect the signal is hidden. I then do a PCA to find wich of the 512 sensors picks up any systematic variation. Then I search for the best fit sinuswave in each of the top sensors (say 5 out of 512).

So main question, how to determine which of the PCA's is the one picking up the true signal without me knowing excactly what the frequency of the signal is? Second question, does the above steps seem reasonable? I am no expert in this, but experiments seems to indicate that with high STN objects of measurements (say a pendulum) it is clearly visible which PCA fits best (amplitude, hz, residuals of fit), but with low STN (say heartbeat), it is not so clear.

Sorry for lengthy text. Thanks alot for any answers!

edit: Running spectral analysis (using different Methods) gives different results depending on wich PC i am analysing. Might there be a spectral Method that does multivariate samples?

edit 2: So while waiting for answers I have chosen to filter out principle components based on the correlation coefficient from a sinus wave fit (r>0.5) and % of data explained (p>0.15). So in other words, the signal i am looking for has to be present in at least 15 out of the 512 sensors (if I understand PCA correctly), and the best fit sinewave has to explain at least 25% of the variation. This works well with test setups in a noisy environment. Question: Is my approach sensible? Should or shouldnt I forego FFT and spectral estimation? Besides visualization, does frequency estimation based on the FFT gain anything over a iterative sinewave fitting, given that I know the signal is sinusoidal?

Thanks alot


texts: Large spike at the center frequency when using Ettus X310

texts:
I am using the Ettus X310 SDR and GNU radio. On the GNU radio schema I have a standard UHD USRP source block connected to a GUI FFT sink - basically I'm just displaying the spectrum of the raw signal from the device.

Parameters for the USRP source block:

sample rate = 1.024 MS/s,
center freq = 101 MHz,
gain = 42 dB,
bandwidth = 130 MHz

The problem that I have is that I always have a big spike (10-15 dB) directly on the center frequency (no matter what frequency I set). I am relatively new to all this so I would appreciate any pointers on how to get rid of the spike.


texts:
Ladies, Gentlemen,
Please let me again employee you with my problems. In following program there is only one sinusoid with frequency 256 Hz amplitude 100,000 and phase 0.67 pi rad. Testing frequencies grid is of 1 Hz step, from 0 Hz to 1000 Hz. Output has maximum amplitude 99,221 and phase 0.670561 pi rad at 256 Hz. I regard it okay given errors in accuracy.

But there are ups and downs in pseudo-frequencies amplitude while my understanding about DFT is that they should going upward till 256 Hz and then downward. Why these ups and downs? I remain with regards.

// Grid is created by variables lowTestFreq and testFreqStep
// both of double type, and integer type M.
// In comments are units of variables.

#include <iostream>
#include "math.h"

int main()
{
    //Unchanged variables.
    double pi = 3.141592653589793; //pure number.
    int n = 0; // sample.
    double xn = 0; //pure number.
    double Re = 0; //pure number.
    double Im = 0; //pure number.
    double xn_cos = 0; //pure number.
    double xn_sin = 0; //pure number.
    double testAmpl = 0; // pure number.
    double testPhase = 0; // rad and then pirad.

    //Sampling frequency and samples of segment.
    int samplFreq = 50000; // sample/second.
    int N = 3000; //sample.

    // Existing sinusoid.
    double ampl_1 = 100000; //pure number.
    double existFreq_1 = 256; // Hz.
    double phase_1 = 0.67 * pi; // rad.
    double unitAngle_1 = (existFreq_1 / samplFreq) * (2 * pi); // rad/sample.

    // Testing grid.
    double testFreq = 0; // Hz.
    double testUnitAngle = 0; // rad.
    double testAngle = 0; //
    double lowTestFreq = 0; //Hz.
    double testFreqStep = 1; // Hz.
    int M = 1000; // pure number.
    int m = 0; // pure number.

    for (m = 0; m <= M; ++m)
    {
        testFreq = lowTestFreq + (m * testFreqStep); // Hz.
        testUnitAngle = (testFreq / samplFreq) * (2 * pi); // rad/sample.

        Re = 0; Im = 0;

        for (n = 1; n <= N; ++n)
        {
            xn = ampl_1 * sin(n * unitAngle_1 + phase_1);

            testAngle = n * testUnitAngle;
            xn_cos = xn *  cos(testAngle);
            xn_sin = xn * -sin(testAngle);

            Re += xn_cos;
            Im += xn_sin;
        }
        std::cout << "Test Freq " << testFreq; // Hz.
        testAmpl = sqrt(Re * Re + Im * Im) / N * 2; //pure number.
        std::cout << "       Test Freq's Ampl: " << testAmpl;
        testPhase = atan2 (Im, Re); // rad.
        testPhase = testPhase / pi; // pirad.
        testPhase = testPhase + 0.5; // pirad. Corrected phase.
        std::cout << ",    and its phase: " << testPhase
            << " pirad." << std::endl;
    }
    return 0;
}


texts: Assessing relative wireless channel quality using Capacity and Condition Number CDFs

texts:
I have Channel capacity analysis figures for two wireless channels show below:

My question:

While ergodic capacity in the first subplots are very clear in conveying the channel capacity, How to make an assessment of channel capacity based on P(Capacity < x) in second subplots and P(Cond. Num < x) in the third subplots for each of the above channels, individually and relatively?


texts: Estimating the PSD (Periodogram): How do you find the Expected Value

texts:
The periodogram $\widehat{P_x}(e^{jω})$ is an estimate of PSD $P_x(e^{jω})$ and for mean-square convergence -

$\lim_{N \rightarrow \infty}E\left\{\widehat{P_x}(e^{jω})\right\}=P_x(e^{jω})$

My question:

If I calculate the expected value, which is according to definition an average, this way

$E\left\{\widehat{P_x}(e^{jω})\right\} = \frac{1}{2π}\int_{<2π>} \widehat{P_x}(e^{jω})dω$

Then, $E\left\{\widehat{P_x}(e^{jω})\right\}$ should have a single value. If this was true, then it would mean that $P_x(e^{jω})$ has the same value for any $ω$ and this clearly is not true. So, where am I mistaking?


texts: I have been studying about the box window and am unable to understand one thing, do we first sample a signal and then window it resulting in up-scaling or downscaling, or is windowing a practical case of sampling itself which is performed on a continuous time signal?


texts: Filter before or after multiplication of two signals?

texts:
I have two signals and want to multiply them, e.g. voltage and current are multiplied to get power. The result (e.g. power) shall be a filtered result, e.g. to see the average.

Is it better to filter (digitally) each signal separately before the multiplication or  to filter after the multiplication?

Intuitively, I would have said that filtering before the multiplication is better for three reasons:

Less numerical problems
Noise of both signals is added by the multiplication, thus filtering before multiplication improves SNR
Multiplication generates modulation products and if you filter out a frequency you do not want to see anyway it also does not cause modulation products.

Am I on the right track here or is there something I am overlooking?

Related: SNR After Multiplying Two Noisy Signals

Edit 1: To answer the question of typical signal spectra: Think of a laboratory grade measurement device like a multimeter. It acquires two signals which can be any signal shape (in time domain) but are band-limited by the input circuitry to avoid anti-aliasing. Noise is typically white noise with a bit of 1/f noise added. SNR is usually 10's of dB.


texts:
I want to calculate averaged power of a time-domain signal by means of its spectrum. I guess Parseval is the right tool.

So I sample a sinus of 100 Hz 10000x within one second.

Unfortunately the sum of the squared samples euqals not the sum of the FFT amplitudes (weighted by the number of FFT bins). Where is the mistake?

# Some python code

import numpy
import matplotlib.pyplot as plt

# Create Time Domain Signal for 1 sec
fs = n = 10000    # Samplingfrequency
ti = numpy.linspace(0,1,num=fs)
sx = 1*numpy.sin(2*numpy.pi*100*ti)

# Calculate spectrum via FFT and account for scaling n/2
# taking the real fft (rfft) only the positive frequencies are calculated
fx = numpy.fft.rfft(sx)/(n/2)
no_of_points = fx.shape[0]

# Calculate RMS for time domains signal + spectrum
parseval_sx = numpy.sum(sx**2)
parseval_fx = numpy.sum(numpy.abs(fx)**2)/no_of_points

print parseval_sx, " equals not ", parseval_fx

Output:

4999.5  equals not  0.000199940012002


texts: How to prepare and plot unequally spaced, irregular data to a contour plot or similar with MATLAB

texts:
I've got a data set of hot-wire measurement velocity amplitudes at a given frequency bin (time data that has already been transformed to the frequency domain and I am just considering data for a given frequency bin) that have been sequentially recorded on specific locations in the xz-plane with a traverse-system. The problem is that the data does not represent a full uniform grid, but is sparsly sampled. Actually, the traversing process was performed in a way to get several locations within a given, more or less rectangular area, while avoiding hitting a circular obstacle within this area. Additionally, the regions close to the obstacle were sampled more densely as compared to regions farther away from the obstacle.

As a result the data set to start my problem with looks like the following  matrix (3 columns and e.g. 400 rows):

x z values (3 columns)

x1 z1 value1
x2 z2 value2
x3 z3 value3
.   .    .
.   .    .
.   .    .

(each value representing the velocity amplitude, e.g. at 500 Hz)

But the problem is:

The data was collected in some wierd order, e.g. the first 6 coordinates (x-component) are [-5; -5; -3; 10; 10; 0.8; ...] and the first 6 coordinates (z-component) are [0.2; 0.3; 5; 3; 2; -1; ...] ...
The data does not represent a uniform and equally spaced grid: all data (for example) lies within the range [-5 <= x <= 10] and [-2 <= z <= 6], but it is possible that there are 8 z-positions for one x-position, while at another x-position there are only 3 z-positions and at several other x-positions there are no z-positions at all (as within the obstacle region).
There is a greatest common divisor with respect to the grid resolution, e.g.    0.1   so all x and z-coordinates are positive or negative multiples of 0.1. This means it would be possible to generate a (probably) huge matrix (xz-plane) that includes all possible locations. I think this could help somehow.

What I basically want:

A simple contourf plot of the data within the xz-plane and areas that have not been recorded (obstacle region) are either interpolated or even better filled with NaNs (or Zeros) or something like that. Actually, I know the coordinates of the obstacle and might just draw it as overlay later.

Do you know how to get the desired result in a convenient manner?

I've been trying for ages now, but I just don't find a solution that solves my problem. I am pretty sure that it can be solved with some combination of reshaping, ndgrid/mesh, sorting or gridded interpolation with built-in (image  processing?) functions but I just don't get it to work.


texts: Understanding DFT-ODD operation in Gardner Efficient Convolution paper

texts:
As a newbie to DSP I am trying to understand and implement an efficient convolution engine as per JAES V43.3 1995/03  - Gardner - Efficient Convolution without I/O Delay.

I'm pretty much there with a working proof of concept, but I am struggling to get my head around appendix two - which talks about DFT-ODD operations and computing spectra using the spectra of two existing blocks of results.

Can someone point me to some additional (readable for newbie) documentations that perhaps covers this topic. This is the last piece in my jigsaw puzzle.


texts: Lowpass then Inverse filtering in python

texts:
I wrote a simple lowpass filter in python to run against lena. Now I'd like to add Gaussian noise to the lowpass filtered data and then run an inverse filter against the lowpass and try to get the original back (well, as close to original). I'm new to programming in python and not quite sure how to add noise and write the inverse.

import matplotlib.pyplot as plt
import numpy as np
import scipy.misc
from scipy import ndimage
import Image

def plot(data, title):
    plot.i += 1
    plt.subplot(2,2,plot.i)
    plt.imshow(data)
    plt.gray()
    plt.title(title)
plot.i = 0

 # Load the data...
img = scipy.misc.lena()
data = np.array(img, dtype=float)
plot(data, 'Original')

#lowpass filter
n=5
kernel = kernel = np.ones((n,n))

plt.show()


texts:
I suppose that OFDM can't modulate analog signal. The basic mathematical description for OFDM is
      $$ v(t)=\sum_{k=0}^{N-1}X_{k}\exp(j2\pi k t /T),\quad 0<t<T$$
In this equation, $X_k$ is a discrete value.
So OFDM can not modulate an analog signal.


texts: Concepts: Ideal and Non Ideal Channel Resulting in Inter Symbol Interference

texts:
I am following the book: Blind equalization and Identification by Zhi Ding and Ye Li. In Chapter 2, the concept of T spaced equalizers is presented. It is mentioned that the output of the channel is $x[k] = \sum_{i=-\infty}^\infty h[i]s[k-i] + w[k]$ where $h$ is the impulse response, $s$ are the symbols / source input and $w[k]$ is iid noise. Then it is defined that *

  "When the channel is non-deal, its impulse response $h[k]$ has more
  than one nonzero component. Consequently, undesirable signal
  distrotion as the channel output $x[k]$ depends on multiple symbols in
  $\{s[k]\}$. This phenomenon is known as inter-symbol interference
  (ISI)."

*

I do not understand why the impulse response has to be zero when it is ideal channel - no ISI. If $h$ is zero then $h[i].s[k-i] = 0$ where $dot$ is the multiplication operator. Or am I mistaken and in the equation $h$ is convolved with $s$ and not multiplied and essentially impulse response will be zero for ideal channel. What is the concept and reason? I do know that convolution in time domain is equivalent to multiplication in frequency domain. But the above representation is in time domain. A brief explanation will work and then I can follow on from the key points mentioned in the answer. Thank you.


texts: Estimating Plane Pose Without Knowledge of Intrinsic Camera Parameters

texts:
Assume a camera has no skew, and has square pixels, that is, its camera-calibration matrix, $K$, looks like:
$$
K = \left[\begin{matrix} \alpha&0&u_x\\0&\alpha&u_y\\0&0&1 \end{matrix}\right]
$$

Intuitively, it seems that it should be possible to recover the pose of a rectangular planar object, $P$, (for example a flat piece of paper) by only knowing the projected image coordinates $\vec{p_i}$ for $i={1, 2, 3, 4}$ of the object's four corners without knowledge of $K$.  To be more precise, I only care about the orientation of the unit normal vector, $\vec{n}$, of $P$, not its location.

Following the answer given in response to this question:  Step by Step Camera Pose Estimation for Visual Tracking and Planar Markers,

We choose world-coordinates such that $P$ is the plane $Z=0$.  We find the homography, $H$, mapping the four-corners of $P$ in world-frame coordinates to its projected coordinates $\vec{p_i}$ for $i=1,2,3,4$.

Then, $H = K[R | \vec{t}]$.  Where $R \in SO(3)$ is the camera's rotation and $t \in \mathbb{R}^3$ is the camera's translation.  If we make the assumption that $K$ is identity (I'll return to this later), then
$H = [\vec{R_1}|\vec{R_2}|\vec{t}]$.  Where $\vec{R_1}, \vec{R_2}$ are the first and second columns of $R$.  The third column of $R$,  $\vec{R_3}$, is equal to $\vec{R_1}\times\vec{R_2}$.  Finally, $\vec{n} = R\left[\begin{matrix}0\\0\\1\end{matrix}\right] = \vec{R_1}\times\vec{R_2}$ and we are done.

Now let's examine the effect of $K$ on our answer $\vec{n}$.  My intuition tells me that $K$ shouldn't matter because varying $\alpha$ will ultimately act as a scalar multiple on our answer $\vec{n}$ which is normalized out since I am only concerned with the direction of $\vec{n}$, not its magnitude.  In addition, it seems that varying $u_x, u_y$ should have the effect of translating the locations of the projected corners $\vec{p_i}$, which should affect $\vec{t}$, but not $\vec{n}$.

Let's see if this is true.  If $K \neq I$, then $H = K[\vec{R_1}|\vec{R_2}|\vec{t}]$ and $$\begin{equation}
\vec{n} = K \vec{R_1} \times K \vec{R_2} = det(K)K^{-T}(\vec{R_1}\times\vec{R_2})\end{equation}$$

In particular, if $\alpha \neq 1$,
$$\vec{n} = \alpha(\vec{R_1} \times \vec{R_2})$$
which is our answer from above times a scalar as expected.

On the other hand, if $u_x, u_y \neq 0$ then
$$det(K)K^{-T} = \left[\begin{matrix} 1&0&0\\0&1&0\\-u_x&-u_y&1 \end{matrix}\right]$$

and
$$\vec{n} = \left[\begin{matrix} 1&0&0\\0&1&0\\-u_x&-u_y&1 \end{matrix}\right](\vec{R_1} \times \vec{R_2})$$

In other words, our choice of camera center has a large effect on our calculations.  This is not intuitive to me.  Isn't our choice of camera center, $u_x, u_y$ completely arbitrary (for example $[0, 0]$ vs. $[w/2, h/2]$) or is there something I am missing?


texts: Method of reconstructing a band-limited signal from discrete samples

texts: I understand current cell phones use digital communications. Given that this industry brings in billions of dollars each year, there is much incentive to get the best performance possible. So what method of reconstructing an audio signal from discrete samples is normally used in a good cell phone these days? It seems there must be methods that work better than the ones I read about in text/reference books.


texts:
From this impulse response:

$$h[n]=a^n u[n]\quad\text{with}\quad\lvert a\rvert<1$$

Applying DTFT analysis equation I obtain this frequency response:

$$H(e^{j\omega})=\dfrac{1}{1-ae^{-j\omega}}$$

But when I plot $\mathrm{abs}[H(e^{j\omega})]$ with MATLAB I obtain something that is a-periodic:

What I am doing wrong? This is my code:

W = [-4*pi:0.1:4*pi];
a = 0.9;
H = 1.0./(1.0-e*a.^(-j*W));
plot(W,abs(H));


texts:
I'm creating an HTML that applies a low-pass Windowed Sinc filter with a window size of 8192 to white, pink, or Brownian noise. There are nine radio buttons, eight of which set a cutoff frequency and one which disables the filter.

I've got the noise generator working, with all three "colours" generating noise with no problem.

I've been looking for a way to implement a Windowed Sinc filter using the Blackman window in JavaScript. I've found one site that seemed to be pretty promising,
but on that site, the event listeners are implemented and the processing itself is done (in a script called "dfilter-0.js") which is extremely large and messy. There are so many unnecessary functions and chains of functions that I can never figure out where anything is or what anything does.

My question is, could I get some help cleaning up the "dfilter-0.js" (none of the online tools I've found worked)? Or alternatively, how does one implement a low pass Windowed Sinc filter with the parameters (cutoff frequency, Blackman window size) and connect it to an AudioContext?

I've found an Audacity plugin that performs the filtering that I want, but it uses the Nyquist language. How would I implement that in JavaScript and HTML?


texts: Is there any free software that can recognize and classify sounds?

texts:
i am making a smart home automation software and i would like to know if there is any software that can recognize noises such as speech ,playing music,phone ringing etc
The exact software i am looking for has already been made by mitsubishi but i cant find the source.
A video by mitsubishi explaining what i need
http://www.merl.com/areas/SoundRecognition/classifier-on-pda.mpeg
also the project's site http://www.merl.com/areas/SoundRecognition/


texts: How can faux correlated data be generated for testing or training?

texts:
There are one or two questions here that ask how to assess data correlation. The data tends to be empirical.

Is there some tool or standardised technique for generating correlated data from scratch?  Some way of generating a sequence with an exact amount of known correlation or autocorrelation? It would allow quantitative testing, training and assessment of correlation measurement tools. Or do people just wing it and knock something up Pythonic?


texts: Angle of arrival (AOA) estimation using FMCW radar using MUSIC algorithm

texts:
I am working with FMCW phased array radar with only upchirps. The range doppler matrix is obtained using the two dimensional fast Fourier transform on multiple chirps.
I want to obtain the angular spectrum per range velocity cell, so have taken range doppler bin for every antenna element and now I have a sample vector rdSnapShot (1 X Nant), where Nant is number of antenna elements present.

For this functionality I am using MATLAB's pmusic (MUSIC algorithm), my question is that as MATLAB's function pmusic requires correlation matrix, is it correct to just multiply the sample vector (1 X Nant) with its complex conjugate and obtain correlation matrix of (size Nant X Nant).

N.B. I have just one range doppler bin for each antenna element.

Providing the MATLAB code:

nTargets = 2;
Rxx_outerproduct = rdSnapShot*rdSnapShot';

[powerSpectrumMUSIC,idxFreq] =pmusic(Rxx_outerproduct ,nTargets,'corr');
% Adjusting the frequency for zero centered angle
idxFreq = idxFreq - max(idxFreq)./2;
% Shifting the response to obtain zero centeric response
powerSpectrumMUSIC=(powerSpectrumMUSIC);
% Converting the frequency scale to direction of arrival scale
idxAngle = asind(idxFreq/pi);
angResponseMUSIC = (db(powerSpectrumMUSIC));


texts: Can an infinite bandwidth signal be an infinite-time signal?

texts:
While decomposing the impulse train using the Fourier series the number of the samples are infinite hence representing a infinite length of train.

But while studying the Nyquist Criteria in the frequency domain it is shown
the expression to have $m$ samples where  $-\infty<m<\infty$

While studying the concept of the Nyquist Criteria in the ISI model where for the system to be interference free the value of the sampling instant must be given as following:

$$h\left(nT_s\right)=\begin{cases}1;&n=0\\0;&n\neq0\end{cases}$$

For more reference here. I cannot understand that the signal is having infinite duration and infnite bandwidth accordingly is it even possible.


texts: How to combine/append mfcc features with rmse and fft using librosa in python 2.7

texts:
I have an audio file say myfile.wav

I trained a neural network based on fft features, and it is giving pretty good results for detecting particular classes of sounds.

I want to expand above experiment to include more sophisthicated features like MFCC along with simpler features like RMSEnergy and so on.

Using librosa,

I can compute mfcc as follows:

import librosa
y, sr = librosa.load('myfile.wav')
mfcc=librosa.feature.mfcc(y=y, sr=sr)
print mfcc.shape
(20,5911)

Similarly  I can compute rmse as follows:

myrmse = librosa.feature.mfcc(y=y)
print myrmse.shape
(1,5911)

Q1.
Is it valid if I transpose mfcc and the myrsme and combine the two using numpy.hstack ?

say

t = mfcc.transpose()
t2 = myrmse.transpose()

So that I now have two arrays : one of size (5911,20) and another of size (5911,1)
and

t3 = numpy.hstack((t,t2))

Q2.
If above operation is valid, how do I  append more features to t3 , Ex:  SNR , FFT  and other audio features.
I'm particularly interested in combining fft with mfcc.

END GOAL
My plan is to have feature data like this :

feat_1  feat_2 ...feat_20  rmse  feat_21 ...feat_n  label
val     val        val      val     val        val      x
val     val        val      val     val        val      x
val     val        val      val     val        val      y
val     val        val

feat_1 .. feat_20 are mfcc related features.
rmse is rms energy from librosa
feat_21 could be some other thing like fft or SNR.

I want to train a ANN and see which combination of features work best for detecting a particular class of sound. And conclude if mfcc's / fft's alone is sufficient for achieving this.

If this is not the right forum kindly direct me to the correct place.

Thank you.


texts: Images similarity measure using j-divergence

texts:
I need to reproduce the result of Xu, Xiaocong, et al, 2016 paper . Mutual information code already was done and work fine. can anyone help me to write Jeffrey’s divergence code? because the I use Jeffrey’s divergence fun in file exchange code but I think it has an issue.

In the aforementioned paper the j-div is better because it has wider distance between correct result and wrong one, but in my case, the opposite MI is better (see the image below)

(up paper result, down my result)

Thanks


texts:
Let just say I understand the first process that is white noise,

  Process 1: If $x(n)$ is a Gaussian random variable and a process
  is formed from a sequence of $x(n)$ and all random variables are
  uncorrelated then this process is white noise.

  Process 2: Now if $x(n) = \alpha $ where now $\alpha$ is Gaussian
  random variable then each process in ensamble is equal to constant for
  all $n$

But I am not able to get what difference does it makes when we form Process 2 it seems same to me. Can any one explain the difference? Also how the ensamble became constant in the second case?


texts:
The null-to-null bandwidth of the main lobe of a BPSK signal is twice the symbol rate.

BPSK transmits one bit per symbol, so the signal bandwidth is twice the bit rate.

That has been shown on a few stack overflow questions. This implies the spectral efficiency is 0.5 bits/second/Hz.

However, a quick Google for the spectral efficiency of BPSK gives an answer of 1 bits/second/Hz. Why?


texts: Obtaining databases for speech emotion recognition analysis

texts: For a research work on speech emotion recognition, I need the database of CASIA and EESDB for simulations. I could not find those; how can I get them? Thanks.


texts:
What is the actual frecuency bandwidth of human speech.
I know there is variation for gender, age and so, but what is the general range... I seem to find different ranges.

Like here http://recognize-speech.com/preprocessing

Is the range defined between 100 hz 8khz

Telephones are made to work in the ranges 300 Hz to 3400 Hz, according to Wiki: https://en.wikipedia.org/wiki/Voice_frequency

The hearing is in range


texts: Difference between $H(\omega)$ and $H(j\omega)$

texts: What is the difference between $H(\omega)$ and $H(j\omega)$, in terms of impulse response. And if so, when do I have to use one or another.


texts: Is there a way to derive an FIR filter using an IIR filter?

texts: I know there's a thread similar to this one here, but the OP is asking the reverse of what I'm trying to find here. I've done some research on the web with very few sources coming up with actual solutions to this problem. What techniques are used to give an approximation of an FIR filter given one or more IIR filters with say the same order?


texts: Does the Fourier series coefficient of AC components remains same if DC component is subtracted form the given signal?

texts:
Suppose a signal is defined by
$
x(t)= \begin{cases}
t & 0\leq t \leq 1 \\
2-t & 1\leq t\leq 2 \\
\end{cases}
$

Since $x(t)$ has even symmetry, I can calculate fourier coefficient as
$$
a_n = \frac{4}{T} \int_0^1 x(t).\cos{n\pi t}.{dx}
$$
I have calculated $$\begin{equation}
a_n = 2\big[\frac{\cos{n\pi} - 1}{n^{2}{\pi}^{2}}\big]\tag{1}
\end{equation}$$

The DC value of $x(t)$ i.e $a_0 = 0.5$. If we subtract DC value we get,

From this we can see that given signal has hidden half wave symmetry in addition to Even symmetry. So we can find fourier coefficient as

$$
a^{'}_n=\frac{8}{T}\int_{0}^{\frac{1}{2}}(t-\frac{1}{2})\cos{n\pi t}.dt
$$

I have calculated $$\begin{equation}
a^{'}_n = 4\big[\frac{\cos{\frac{n\pi}{2}} - 1}{n^{2}{\pi}^{2}}\big]\tag{2}
\end{equation}$$

My question is, shouldn't $a_n$ and $a^{'}_n$ be equal for $n\neq0$ ?


texts: Demodulating wideband FM signal in Octave from recorded baseband IQ samples

texts:
Assuming I have a baseband recording of 8 bit IQ samples at 2.4MSps centered on an FM station is it correct to say that after low pass-filtering the complex signal to 240KHz, taking discrete derivative of phase and computing the envelope of the result this envelope should be the baseband audio signal?

Will this approach work for FM stereo?

I'm trying to confirm that my understanding is correct in Octave based on IQ samples recorded from RTL-SDR tuned to an FM station, but I'm getting white noise.

pkg load signal;
% 2.4MSps IQ baseband signal of an FM station from RTL-SDR recorded in SDR#
samples = audioread("C:\\Users\\remote\\Downloads\\sdrsharp-x86\\SDRSharp_20191228_175834Z_101200000Hz_IQ.wav", "native");
csamples = arrayfun(@complex,samples(:,1),samples(:,2));
% low pass filter and decimate to get 480KSps (240KHz Widband FM)
csamples_filtered = decimate(csamples, 5);
% descrete derivatives of phase
filtered_derivative = conv(angle(csamples_filtered), [-1,1]);
% center amplitude around 0
filtered_derivative = filtered_derivative - (max(filtered_derivative)+min(filtered_derivative));
% rectify, low pass filter and decimate to get the envelope (48 KSps)
envelope = decimate(sqrt(filtered_derivative.*filtered_derivative), 10);
% scale to 0..255
sound_normalized = sound_normalized * (255/max(envelope));

player = audioplayer(sound_normalized, 48000, 8);
play(player);

The FFT plot of the original data in csamples (centered on Fc) corresponds to spectrum that is shown in SDR Sharp I recorded it in (and listened to the station), so I think the input data is good.


texts: What Does Multi Scale Image Analysis Mean?

texts:
I'm reading an article about multi-scale representation of image, and it is said that, convolving the image with a Gaussian kernel at different $\sigma$, then different scale representation is created.

This really confuses me, because I think the image is just blurred by Gaussian with different $\sigma$, each blurred version still has the same image size, so I don't think they're scaled, unless each blurred is downsampled to a smaller image size, which means I think different size representation is truely the


texts: Sinc interpolation of pure sine wave sampled just at Nyquist frequency

texts:
Following this question: Shannon-Nyquist theorem reconstruct 1Hz sine wave from 2 samples

could you explain the algorithm to apply for sinc interpolation to avoid the "sawtooth" effect due to linear interpolation?

(It seems to me that the https://en.wikipedia.org/wiki/Whittaker%E2%80%93Shannon_interpolation_formula  Shannon-Whittaker formula would be the one suitable for this?)


texts: MATLAB Implementation of Karplus Strong algorithm with filter function?

texts:
I want to implement following function: y = ksalgrithm(x, alpha, M, Nout) where

x is the input vector with length Nout
alpha=0.99
Nout is the length of the output vector.

The algorithm should do the convolution of $y[n]=(h_1\star x)[n]$, where $h_1[n]$ is the inverse of the $z$-transformation of $H_1(z)$. The input signal should have values $l = T_0 \cdot f$ and filled up with zeros until the vector has the length of Nout.

I want to use MATLAB's filter.

Nout = length(x);
M = T_0*fs;
alpha = 0.99;
function [ y ] = ksalgorithm( x, alpha, M, Nout )
    syms z;
    H = filter(1,1-alpha*z^(-M),x);
    H_1 = 1/(1-alpha*z^(-M));
    h_1 =  iztrans(H_1);

    y = zeros(x);
    ...

I get always this error?

Undefined function 'filter' for input arguments of type 'sym'.
Error in ksalgorithm (line 23)
    H = filter(1,1-alpha*z^(-M),x);


texts: Intuitive interpretation of Laplace transform

texts:
So I am getting to grasps with Fourier transforms. Intuitively now I definately understand what it does and will soon follow some classes on the mathematics (so the actual subject). But then I go on reading about the laplace transform and there I kind of lose it. What is the moment of a signal? Why is the fourier transform a special case of the laplace transform? How can I come to grips with the Laplace transform?

Ive looked at these sources before I asked this question:

What is meant by a system's "impulse response" and "frequency response?"

How to distinguish between the different frequency domains?

Amplitude vs Frequency Response

Why is the Fourier transform so important?

http://en.wikipedia.org/wiki/Laplace_transform


texts:
I am implementing MUSIC algorithm in Verilog, and I need to implement the blocks in following order
$$\text{Input signal}{\longrightarrow}\text{Co-variance/Correlation matrix}{\longrightarrow}\text{EVD}{\longrightarrow}\text{Peak search}$$
But I am confused because some papers refer the second block as correlation and some it as co-variance.


texts: what is the reason for the Change in DFT analysis after downsampling

texts:
firstly the original signal's dft analysis is this the signal is sammpled at 1000Hz

the signal is then downsampled under different rates such as 25Hz

this is the dft analysis of the signal downsampled under 25Hz

what is the reason for this phenomenon?


texts: How to choose order and cut-off frequency for low-pass Butterworth filter?

texts:
When using IMU data (accelerations and angular velocities), many papers use a 4th order low-pass butterworth filter with a cutoff frequency of 5 Hz. However, none of them really justify why those filters have been chosen.

Can anyone please explain how the order and cutoff frequency are chosen for a filter of this nature? My data has a sampling frequency of 125 Hz


texts: GCC-PHAT (Generalized cross correlation MATLAB)

texts:
I've been trying to implement generalized cross correlation with a PHAT weighting function for a while now, and cannot get it to work.  I have tried correlating using MATLAB's xcorr.m file, and it works for getting a correct delay tau (on simulated sinusoidal signals).

The code for the signals are thus:

Fs = 8000;
dt = 1/Fs;%0.125e-3
f1 = 100;
tdelay = 0.625e-03;%try different values
t3 = (0:dt:(1)-dt)';
x3 = cos(2*pi*f1*t3);
x4 = cos(2*pi*f1*(t3-tdelay));

As can be seen in part of the source code for xcorr.m, the cross correlation is implemented thus:

%Transform both vectors

X = fft(x,2^nextpow2(2*M-1));
Y = fft(y,2^nextpow2(2*M-1));

% Compute cross-correlation

c = ifft(X.*conj(Y));

According to definitions of GCC-PHAT, the only addition I needed to make was to divide the product by its own magnitude, before taking the ifft.  Here is my version with this change.

%Transform both vectors
X = fft(x,2^nextpow2(2*M-1));
Y = fft(y,2^nextpow2(2*M-1));

% Compute cross-correlation

R = X.*conj(Y);
c = ifft(R./abs(R));

However I always end up with a tau of zero with the PHAT weighting!  On closer examination of the array produced as a result of this division, it seems as if the first value of R is a real value (with no imaginary component), and so when divided by its magnitude it becomes 1.  All the other values in array R are complex, so do not end up as 1 when divided by their own magnitude and so end up with a value of <1.

This can be seen below, for the 1st 10 values of R.

    K>> R(1:10,1)

ans =

  0.000000000000000 + 0.000000000000000i
 -0.494299608718696 - 0.003002230689022i
 -0.002678647083223 - 0.000032538742345i
 -0.488954228290329 - 0.008909374553649i
 -0.010656518992354 - 0.000258902698589i
 -0.478379290671260 - 0.014528074329782i
 -0.023760667475633 - 0.000865926459320i
 -0.462803929640386 - 0.019677623220519i
 -0.041707017319469 - 0.002026674993917i
 -0.442565618721743 - 0.024194329448597i

K>> abs(R(1:10,1))

ans =

   0.000000000000000
   0.494308725968464
   0.002678844707371
   0.489035391682370
   0.010659663580139
   0.478599844010494
   0.023776441018801
   0.463222070011989
   0.041756229537848
   0.443226457301486

K>> R(1:10,1)./abs(R(1:10,1))

ans =

  1.000000000000000 + 0.000000000000000i
 -0.999981555555690 - 0.006073594357736i
 -0.999926227844417 - 0.012146557900713i
 -0.999834033705084 - 0.018218261306199i
 -0.999705001216859 - 0.024288074069393i
 -0.999539169638280 - 0.030355367874845i
 -0.999336589392997 - 0.036419515378054i
 -0.999097322000240 - 0.042479891383435i
 -0.998821440083941 - 0.048535871565711i
 -0.998509027227826 - 0.054586834901284i

As can be seen above, the largest value ends up being at the 1st index, when we divide R by its own magnitude.  So once the ifft is taken, the highest value is ALWAYS at the beginning of the array, which gives a lag and timediff of zero...even when I set the delay between the 2 identical signals at the beginning to a nonzero value (e.g. delay = 0.75e-03).

What am I doing wrong?? Any help appreciated.

Rory

P.S. if anyone is wondering why I'm bothering with the PHAT weighting, its because it should give far better results (in theory) in a real life scenario, for TDOA.


texts: Question about symbol detection in ISI channels (use of APP decoder in Matlab)

texts:
I am studying the problem of digital symbols detection in ISI channel. In the figure below, the output symbol is represented as a weighted sum of input symbols.

I am trying to use MAP decoder comm.APPDecoder to estimate the transmit symbols. To construct the trellis of such figure, I am using the Matlab function poly2trellis. The function can be used as follows:

trellis = poly2trellis(ConstraintLength,CodeGenerator)

My question is how can I enter (or consider) the real-valued coefficients $\Gamma_{-L_2}, ..., \Gamma_{L_1}$ in the second input argument CodeGenerator. If these coefficients are binary, then it is straightforward to convert this binary sequence to octal. But, I have no clue with real-valued coefficients. Any hints?


texts: Is there any relationship between the FFT and linear prediction?

texts: Is there any relationship between the FFT and linear predictive methods?  Can an FFT result or the input to an FFT be modified to do non-circular prediction/extrapolation from the FFT results.


texts: What type of correlation is this equation?

texts:
In a book [1], gives an equation (equation 5.31 in the book) for estimating frequency offset in OFDM systems proposed by Classen. the procedure is:

  first, two OFDM symbols, $y_l [n]$ and $y_{l+D} [n]$ , are saved in a memory after synchronization. then, the signals are transformed into ${Y_l [k]}_{k=0}^{N-1}$ and ${Y_{l+D} [k]}_{k=0}^{N-1}$ via FFT, from which pilot tones are extracted. After estimating CFO from pilot tones in the frequency domain, the signal is compensated in the time domain.

Here is the equation:

$ \hat \epsilon_{acq} = \dfrac{1}{2\pi \cdot T_{sub}} \underset{\epsilon}{max} \left \{ \left \vert \sum_{j=0}^{L-1} Y_{l+D}[p[j],\epsilon]\cdot Y_{l}^{*}[p[j],\epsilon]\cdot X_{l+D}^{*}[p[j]]\cdot X_l[p[j]] \right \vert \right \}  $

where $L$, $p[j]$, and $X_{l}[p[j]]$ are the number of pilot tones, the location of the jth pilot tone, and the pilot tone located at $p[j]$ in the frequency domain at the $l^{th}$ symbol period, respectively.

From the equation, It's supposed to find the normalized CFO, $\epsilon$. but, the parameter we are finding is involved in the equation $[P[j],ε]$.
My question is: how this parameter we're supposed to estimate is in the middle of the equation? What type of correlation is this? Or did I misunderstand the equation? Please explain.

[1] Yong Soo Cho et al., "MIMO-OFDM wireless communications with MATLAB", John Wiley & Sons (asia), 2010


texts: MFSK vs FSK: Why BER decreases increasing the symbols

texts:
My question is related to the performance of 2-FSK modulation scheme against M-FSK (M > 2).
From the SNR vs BER curve it seems that as the number of tones (frequencies) increases the BER decrease for the same SNR compared to 2-FSK(pics attached). (I assume maintaining the same symbol duration)

As the number of tones increases the bandwidth increases but i can't figure out why the performance of MFSK also increases.
At first glance i would say that the resilience to the noise should be the same as the 2-FSK.

Which is the benefit of using 2-FSK instead of M-FSK (M>2) given a fixed available bandwidth?

Peter


texts: How can I plot frequency vs magnitude of wavelet transform?

texts:
I am running Morlet continuous wavelet transform. I have got wscalogram of signal and now I want to plot freq-magnitude like the following picture., but I don't know how to do it:

I've used scal2freq MATLAB function to convert scales to pseudo frequencies. Also I have some frequencies in my signal that have a large damping ratio(4%) so they are not well visible in plot. How can I exaggerate these highly damped modes?

I am using MATLAB, here is my codes:

% Import the text4.txt to matlab workspace. and save it under name "data"
t=linspace(0,30,301);
Fs=ceil(inv(t(2)-t(1)));
x=data(:,4); % use x=data(:,3),x=data(:,5) too. first column is time,second is refrence
wname = 'morl';
scales = 1:1:256;
coefs = cwt(x,scales,wname,'lvlabs');
freq = scal2frq(scales,wname,1/Fs);
surf(t,freq,abs(coefs));shading('interp');
axis tight; xlabel('Seconds'); ylabel('Pseudo-Frequency (Hz)');
axis([0 30 0 1 0 60])
xlabel('Time'); ylabel('Frequency');
figure;
sc=wscalogram('image',coefs,'scales',freq,'ydata',x);
hold on

My current plot looks like this:


texts: Phase spectrum of the DTFT of a rect pulse using FFT

texts:
I'm trying to use the FFT to plot the phase spectrum of the DTFT of a sampled rect impulse function.
Here's the Matlab code I've written.

fs=10; % sampling frequency Hz
Ts=1/fs;
T=4; % length of time vector in seconds
t=-T/2:Ts:T/2-Ts; %vector of 40 time samples -2:0.1:1.9
y=rectpuls(t); %unit length rect impulse
z=fft(y,10000).*Ts;
%fft multiplied by Ts to get rid of the 1/Ts factor in the spectral replicas
%of the DTFT of the signal
f=linspace(0, (length(z)-1)/length(z),length(z));
figure
plot(f,abs(z));
grid on

%plot of the phase of z as it is, so with the time sequence
%considered to be shifted by length(y)/2 samples by fft.
figure
plot(f,unwrap(angle(z)),'r');
grid on

%I try to get rid of the linear phase delay component to get the phase
%spectrum of the DTFT of the signal in its frequency period.
%I divide each fft sample by its own exponential delay factor.
delay=length(y)/(2*length(z));
expd=zeros(1,length(z));
for index=1:1:length(z)
expd(index)=exp(-1i*2*pi*(index-1)*delay);
end
z2=z./expd;
figure
plot(f,unwrap(angle(z2))./pi,'r');
grid on

Now the problem is that this phase plot is wrong, the DTFT is the sum of real sinc functions so its phase spectrum con only assume k*pi values.
There's still a linear component in the phase plot and the slope is positive so the exponentials that I'm using are over compensanting the delay.
I have found out that the problem is that the fft function introduces a linear phase delay as if the sequence was shifted by 19.5 samples instead of 20 samples as it should be in my example.
Infact if I set in my code delay=19.5/10000 I get the correct phase plot (or at least one that makes sense :D):

So what's going on here, is it a precision problem?
Are there mistakes in my Matlab code and/or my dsp theory formulas?
Thank you!!

Unfortunately I can't link more than two images in the post :/


texts: how to detect short pause in a speech with noise?

texts:
I would like to count the number of short pauses in a speech and know the duration of the pauses by loading the signal into Matlab but I am very new to Matlab.
Can anyone help and give me some ideas about this?


texts: What do you call the random Gaussian vectors in compressed sensing?

texts:
Let $(e_i)_{1\leq i\leq M}$ be vectors with zero-mean i.i.d. entries from a Gaussian distribution. In a compressed-sensing setup, I have observed the collection of scalar products $(\langle x\, |\, e_i\rangle)_{1\leq i\leq M}$ for an unknown vector $x$.

In the compressed sensing literature, what is the name for the collection of vectors   $(e_i)_{1\leq i\leq M}$?


texts: Efficient algorithm for computing area under convolution?

texts:
Consider $N$ discrete signals $x_1(n), x_2(n), \ldots,x_N(n)$ each a bounded support of size $M$.

To convolve them, we can zero-pad each of them, multiply their FFTs, and take the result's inverse FFT.
I believe this requires taking $N$ FFTs of length $N(M-1)+1$ each, taking $O(MN\,\log(MN))$ time.

However, what I need isn't the full convolution. I only need the area under the convolution above zero:
$$\sum_{k=0}^\infty y(n)$$
assuming $y$ is the convolution of our signals:
$$y = x_1 * x_2 * \ldots * x_N$$

Is there any algorithm that allows me to compute this value with a better time complexity than that required for computing the full convolution, or are the problems equally difficult?


texts: kernels to compute second order derivative of digital image

texts:
For an image $I$, its first order derivatives can be computed using several oprators, such as
$$K_{sobel} = \left[ \begin{array}{ccc} -1 &0 &1 \\ -2 &0 &2 \\ -1 &0 &1 \end{array}\right]$$
$$I_{x} = I * K_{sobel},\ I_{y} = I * K_{sobel}.T$$.
And to compute the second order derivatives $(I_{xx}, I_{yy}, I_{xy})$, I can apply $K_{sobel}$ twice on $I$, however, I wonder are there any kernels which can be applied once to get the second order derivatives.

Here is what I tried, I made up 2 kernels,
$$K_{xx} = \left[ \begin{array}{ccc} 0 &0 &0 \\ 2 &-4 &2 \\ 0 &0 &0 \end{array}\right]$$
$$K_{yy} = \left[ \begin{array}{ccc} 0 &2 &0 \\ 0 &-4 &0 \\ 0 &2 &0 \end{array}\right]$$
and $I_{xx} = I * K_{xx},\ I_{yy} = I * K_{yy}$, they look pretty much the same as twice-applied-sobel results, but I failed to find a good kernel for $I_{xy}$, I need your help, any idea is appreciated.


texts:
I was reading a paper about Gist features. I was interested about its code implementation which i found it here). Into LMgist.m there is a function called createGabor, this is the code:

function G = createGabor(or, n)
%
% G = createGabor(numberOfOrientationsPerScale, n);
%
% Precomputes filter transfer functions. All computations are done on the
% Fourier domain.
%
% If you call this function without output arguments it will show the
% tiling of the Fourier domain.
%
% Input
%     numberOfOrientationsPerScale = vector that contains the number of
%                                orientations at each scale (from HF to BF)
%     n = imagesize = [nrows ncols]
%
% output
%     G = transfer functions for a jet of gabor filters

Nscales = length(or);
Nfilters = sum(or);

if length(n) == 1
    n = [n(1) n(1)];
end

l=0;
for i=1:Nscales
    for j=1:or(i)
        l=l+1;
        param(l,:)=[.35 .3/(1.85^(i-1)) 16*or(i)^2/32^2 pi/(or(i))*(j-1)];
    end
end

% Frequencies:
%[fx, fy] = meshgrid(-n/2:n/2-1);
[fx, fy] = meshgrid(-n(2)/2:n(2)/2-1, -n(1)/2:n(1)/2-1);
fr = fftshift(sqrt(fx.^2+fy.^2));
t = fftshift(angle(fx+sqrt(-1)*fy));

% Transfer functions:
G=zeros([n(1) n(2) Nfilters]);
for i=1:Nfilters
    tr=t+param(i,4);
    tr=tr+2*pi*(tr<-pi)-2*pi*(tr>pi);

    G(:,:,i)=exp(-10*param(i,1)*(fr/n(2)/param(i,2)-1).^2-2*param(i,3)*pi*tr.^2);
end

if nargout == 0
    figure
    for i=1:Nfilters
        contour(fx, fy, fftshift(G(:,:,i)),[1 .7 .6],'r');
        hold on
    end
    axis('on')
    axis('equal')
    axis([-n(2)/2 n(2)/2 -n(1)/2 n(1)/2])
    axis('ij')
    xlabel('f_x (cycles per image)')
    ylabel('f_y (cycles per image)')
    grid on
end

If you run this function with the input parameters or = [8 8 8 8] and n = 192 (as in the original implementation) it should create this image:

This image shows the contours of the gabor filters in the frequency domain for z = 0.6, z= 0.7.

Unfortunately this gabor filter implementation is different from what i've seen so far. My questions are:

Look at the matrix param, the second column definitely seems to be
the frequency scale, where fmax = 0.3 and c = 1.85 is the constant.
The fourth column definitely seems to be the orientation scale. The first column is a constant of value 0.35; the only reference i've found in that paper about 0.35 is this sentence (page 8, right column):

  The model of Eq. (7) provides correct fitting for all the eight
  categories for frequencies below 0.35 cycles/pixel (as noise and
  aliasing corrupt higher spatial frequencies, see Fig. 4)

So what the first and third column represent?
Which brings me to this line of code:
G(:,:,i)=exp(-10*param(i,1)*(fr/n(2)/param(i,2)-1).^2-2*param(i,3)* pi *tr.^2);
this implementation is very different from what i've seen so far, i
don't think i recognize the formulation. In particular the second
term of the exponential: 2*param(i,3)* pi *tr.^2 controls the
orientation scale and also the particular "stretch" of the gaussian
in the frequency domain as shown in the picture, which I think is
used to get a better filter overlap in the frequency domain. Can
somebody clear up this implementation?


texts: pitch detection: which is the most robust way to distinguish voiced from unvoiced audio

texts:
I am trying to write a voice/whistle driven synth but I am struggling with pitch detection. As a premise I must say I am quite skilled and experienced at dsp coding plus I have pretty read almost everything I could find in literature about this very thorny topic and I understood that a single robust solution does not exist. I tried and invented every possible and fancy approach in frequency domain based on analysis of the harmonic series or HPS and variations thereof, but I am still plagued by octave errors,  plus the lack of a robust criterion to distinguish pitched from unpitched frames (i.e containing noise). At the end I decided to give up spectral approaches and use time domain autocorrelation (by FFT, squared magnitude plus IFFT) which turned out the most reliable approach of all even if still not 100% error prone. But again, I cannot find any reliable and robust rationale to decide when a frame is voiced/pitched and when it contains noise or transients and therefore the detected highest peak in autocorrelation has to be deemed invalid and discarded.

I found the sourcecode of a plugin called Autotalent which I judged interesting,  the author there relies on a "confidence" value obtained by multiplying thr height of the highest peak in autocorrelation by the value at the corresponding position of the autocorrelation of the window used to window the processed frames but I really failed to understand the logics behind that. I tried for curiosity to use the same trick but what I get is a value which seems completely unrelated to the fact a frame is pitched or unpitched.

I kindly ask where should I start looking for a solution which is not,  possibly,  too complicated (everything is meant for realtime operation). Any hints ?  Do you know any robust method to distinguish pitched and unpitched frames (either in time or frequency domain ) ? Thanks


texts: Can all systems be described by differential/difference equations?

texts:
All LCCODEs (linear, constant-coefficient ordinary differential equation) with initial rest conditions describe causal LTI systems (is this a correct statement?)

Is the converse true, i.e. can all causal LTI systems be described by LCCODEs with initial rest? If not, what are the necessary/sufficient conditions?

It seems trivial that arbitrary differential/difference equations describe some sort of system (is that also true?)

Again, is the converse true, i.e. can we describe any arbitrary systems with (possibly
nonlinear/variable-coefficient/partial, without restriction on initial
conditions) differential equations? If not, what are the necessary/sufficient conditions?


texts: How to prove energy is preserved in sampling using Parseval's Relation?

texts:
I have a continuous signal,xc(t) which is bandlimited to =pi/T ,
and x[n] = xc(nT) are samples of this signal.
I need to show there is a Parseval's relation between the original signal and the sampled version.That is I need to prove that energy will be conserved!
Can anyone please help?


texts: Improve coverage of BackgroundSubtractorMOG2 and findContours

texts:
I'm trying to track players on a soccer field, using a fixed camera. The picture below shows the results of BackgroundSubtractorMOG2, the second picture below the results of findContours on the results of BackgroundSubtractorMOG2.

As you can see there are a lot off small contours detected. How can a I alter the results of BackgroundSubtractorMOG2 so one shape (a player) results in one contour I can track ?

This is the code I used:

import numpy as np
import cv2

cap = cv2.VideoCapture("images/Keeper.mov")
history = 30   # or whatever you want it to be
accelerate = 5

kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(3,3))
fgbg = cv2.BackgroundSubtractorMOG2()

while(1):
    for i in (1, accelerate):
        ret, frame = cap.read()

    fgmask = fgbg.apply(frame, learningRate=1.0/history)

    fgmask = cv2.morphologyEx(fgmask, cv2.MORPH_OPEN, kernel)
    cv2.imshow('frame',fgmask)

    h, w = fgmask.shape[:2]

    contours0, hierarchy = cv2.findContours( fgmask.copy(), cv2.RETR_TREE,
    cv2.CHAIN_APPROX_SIMPLE)
    contours = [cv2.approxPolyDP(cnt, 3, True) for cnt in contours0]

    vis = np.zeros((h, w, 3), np.uint8)
    levels = 17
    cv2.drawContours( vis, contours, (-1, 3)[levels <= 0], (128,255,255),
            -1, cv2.CV_AA, hierarchy, abs(levels) )
    cv2.imshow('contours', vis)

    k = cv2.waitKey(30) & 0xff
    if k == 27:
        break

cap.release()
cv2.destroyAllWindows()


texts:
I'm trying to understand how I can start from the CTFT of a signal and end up with a DTFT.

For example if I take a basic example:

$$\begin{aligned}
x(t) &= \cos(\omega_x \cdot t) = \frac{1}{2} \cdot \left( e^{j\omega_x t} + e^{-j\omega_x t} \right) \\ \implies X(\omega) &= \pi \cdot (\delta(\omega - \omega_x) + \delta(\omega+\omega_x)) \\ x_c(t) &= x(t) \cdot \sum_{n=-\infty}^{\infty}\delta(t - nT_s) = \sum_{n=-\infty}^{\infty}x(nT_s)\delta(t-nT_s) \\ \implies X_c(\omega) &= X(\omega) * \left( \omega_s \sum_{n=-\infty}^{\infty} \delta(\omega - n\omega_s)\right) \\
\end{aligned}$$

$$ X_c(\omega) = \omega_s \pi \sum_{n=-\infty}^{\infty} \left(  \delta(\omega - \omega_x - n\omega_s) +  \delta(\omega + \omega_x - n\omega_s) \right)
\tag{1} \label{1} $$

From there I'm lost and everything crumbles. I'm only trying to get the DTFT of a cosine which is:
$$
\cos(\Omega_0 n) \Longleftrightarrow \pi \sum_{n=-\infty}^{\infty} \left(  \delta(\Omega - \Omega_0 - n2\pi) +  \delta(\Omega + \Omega_0 - n2\pi) \right)
\tag{2}
\label{2}
$$

How can I obtain $\eqref{2}$ starting from $\eqref{1}$?

I hope what I'm trying even makes sense. After all the DTFT with infinite period is the CTFT so I suppose there's a link we can make between these equations?

Thanks


texts: Precompute parts of an FFT on streamed data inputs

texts:
I've implemented Radix-2 & Radix-4, DIT & DIF variants of the cooley-tukey FFT algorithm and they perform well on test data.

Now for a project, my algorithms will be applied to a real-time data stream of complex values. My current approach is window-based, meaning that I'm creating a buffer that stores values in a 2d array for a certain amount of time and after that, the FFT is applied to the matrix stored in the buffer.

As far as my understanding of the algorithm goes (having implemented it and read some papers on DFT/FFT in general but no prior experience in DSP), there is no way of speeding this up by doing some computation live as the data points arrive in the stream connector?

I was wondering if this is correct since I did not find any literature/paper that does this. Any hint is appreciated.

Edit: Example for clarification

Say my input is a Matrix of size N x N of complex numbers and I'm getting the values one point at a time from the data stream, can I do anything more than just buffering them in a 2d array and doing the FFT after a certain time.


texts:
I am writing a lens correctionion code and the process should do several things:

Inverse fisheye projection to obtain rectilinear image
Barrel/pincushion correction
Vignetting correction
Lateral chromatic aberration correction

The problem is - in what order apply these steps?

Should the image be first "de-fished" or the chromatic aberration, vignetting etc. affects already formed fisheye image?

I am not familiar with optics so maybe any ordering is okay as long as correct models/parameters are used - but I would like to go for simplest solution in that case - e.g. if chromatic aberration is most easily corrected in rectilinear images, then de-fishing should go first.


texts: set of signals containing time and amplitude was given to me and asked to process the signal with the help of FFT function. I need help in converting the signal to frequency domain and to extract the feature from the signal so that I can compare the test signal with the given signal so the signal is a fault signal or not


texts:
I hope, you are doing well. I want to ask question about baseband conversion.
Background:

I generated two AM signals in MATLAB, which one of them is shifted in time.
I multiply the signals with $\cos(2\pi f_c t)$ in order to obtain IF signal.
I filtered the IF signals with $f_c$ cut-off frequency.

In this way, I aim to downconverted to signal to baseband.
My problem is, when I try to cross-correlated the baseband signals, my result peaks at zero.
I expect to see the delay, but I could not. Am I missing something?
Notes:
I am exact sure, my cross-correlation function works perfectly. I tried this function before baseband process.
Also, in order to obtain both I and Q components, I try to multiply my signals with sin(2pifc*t). Nothing changes.
MessageSignal: Audio Signal.


texts: How to distinguish signal and noise eigenvalues of autocorrelation matrix for MUSIC algorithm?

texts:
MUSIC (multiple signal classification) algorithm is introduced in Schmidt, Ralph. "Multiple emitter location and signal parameter estimation." IEEE transactions on antennas and propagation 34.3 (1986): 276-280.

An important part of this algorithm is to distinguish eigenvalues froma sample autocorrelation matrix ($R_{N}=\frac{1}{N}\sum{XX^{H}}$) into two groups: noise eigenvalues and signal eigenvalues.

Eigenvalues are plotted as below. A few nontrivial signal eigenvalues and lots of trivial (zero or close to zero) noise eigenvalues:

The figure above is in a paper Quinlan, Angela et al. “Model Order Selection for Short Data: An Exponential Fitting Test (EFT).” EURASIP J. Adv. Sig. Proc. 2007 (2007): n. pag.. The paper is telling how to distinguish signal eigenvalues and noise eigenvalues. The authors use exponential fitting test (EFT), which is not that simple, but seems effective.

On the other hand, I've found two implementations of MUSIC algorithm in Python and MATLAB. Interestingly, they've used very simple approach: calculate a difference (i.e., a slope) between adjacent eigenvalues and find an index of the most largest difference. Then, noise eigenvalues are with indices larger than the index of the largest difference.

Intuitively, EFT can coincide with the largest difference approach. However, it is just an intuition, and does not include an analysis or reasoning.

Is it okay to use the largest difference approach to distinguish signal eigenvalues and noise eigenvalues? If so or if not, what is a good explanation of it?


texts: Second order sections and normalized factor

texts:
I use matlab function tf2sos to convert digital filter transfer function data to second-order sections form:

  [sos,g] = tf2sos(b,a)

Taken from the help menu of the matlab:

  G is a scalar which accounts for the overall gain of the system

My question is: Because G is the overall gain if i draw the block diagram is going to be something like below (i use office visio for the pictures, if you press view picture the resolution is good but in this thumbnails they look awfull, sorry about that) :

But if i want yo design my system like:

how can i find the gain between every node(b1,b2...etc)?


texts:
I'm a SP student and still a beginner in the field.

When we were learning about window functions, my professor said "you're gonna almost always want to window your data but don't window a transient signal." It's a (semi-)online course with lecture videos so I couldn't ask for an explanation.

Can anyone explain what a transient is and why we wouldn't window it?


texts:
I would like to experiment with some input from a microphone and am receiving a 2 channel, 512 samples buffer in real time

I know the signal could be passed through a low pass filter, and then windowed before the FFT.

The sample rate is 44100.00, what low pass filter is needed for this, does the analog to digital conversion not run it through a low pass filter before converting it?

How exactly does a window apply to real time sample input of 512 samples? Is a window not the same as a filter?

I have also read that the signal can be run through the FFT twice, why exactly is this done and is it necessary for this data?

I would like to slowly build this up to do spectral analysis etc. Is this done directly on the FFT bin data?


texts:
How do those matrices work? Do I need to multiple every single pixel? How about the upperleft, upperright, bottomleft and bottomleft pixels where there's no surrounding pixel? And does the matrix work from left to right and from up to bottom or from up to bottom first and then left to right?

Why does this kernel (Edge enhance):

turns into this image:

?


texts: What's the relation between frequency band of $X(j\omega)$ and $\Phi_{xx}(j\omega)$?

texts:
in which:

$x_{c}(t)$ is a continuous-time signal
$X(j\Omega)$ is the Fourier Transform of $x_{c}(t)$
$\Phi_{xx}(j\Omega)$ is the Power Spectrum Density of $x_{c}(t)$ which defined as Fourier Transform of auto-correlation of $x_{c}(t)$.

in other words I want to know when it is said that a signal is band-limited which is the case? band-limited according to $X(j\Omega)$ or $\Phi_{xx}(j\Omega)$?
or if it's a relationship between two cases what it is.
thanks.


texts: What is the optimal way to select by color?

texts:
I need to select all pixels in a image that are similar in color to a given one.

My current approach is to get the hsv values for this pixel and calculate the distance to each pixel in the image being

distance = abs(ref.h-pix.h) + abs(ref.s-pix.s) + abs(ref.v-pix.v);

And select the pixels bellow certain distance.

Is there any other preferred method?


texts: Template matching or object recognition?

texts:
Problem: I have a photo of an object (a manufactured part like the attached photo below), using my Andoird phone camera I want to verify if the object in camera preview matches to the template or not. (in other words, is it the same part as the template or not)

I can make the user to move the camera in order to have similar view of the template in camera preview as the template however there will be different noise level and/or lighting and maybe different background.

Question:  What do you recommend me to use for solving this problem? I was thinking of Canny edge extraction and then matching the camera frames towards the canny edge extract from template? is this a good idea? if yes would you please tell me how can I implement this? any resources? samples? (I can do the Canny edge extraction but couldn't find a way to do the matching)

if Not a good idea then what do you recommend?

Things I have tried:

Feature Extract and Matching: I used few different extractor and matcher implementations from OpenCV and my app is working and drawing the detected feature points and matches, etc. however being a beginner with image processing I cannot make sense of the result and also how to know what is a match. any idea, help, good resources?
Template Matching: I used OpenCV template matching however the performance was horrible and I decided that this cannot be the solution.


texts: Finding period of a square wave with varying sampling frequency

texts:
I have a square wave (0-1.8V) with a varying sampling frequency (from a circuit simulator). It is also not a perfect square wave (the high and low signal could be very close to but not precisely zero or 1 Volt) and the fall time is slightly slower than the rise time.

Another factor is that I will retake the data and the period (and possibly duty cycle as well) will change when I alter another parameter, so I am not working with a fixed number of cycles per window.

I tried autocorrelation but it seems that it gave me double or so the actual period. And fft seems to need a fixed sampling rate.  What would be the best solution to finding the period of this square wave?


texts: Bilateral Laplace transform and existence of Fourier transform

texts:
I was reading from Athanosios Papoulis' "The Fourier integral and its applications." and they referenced the bilateral Laplace transform and Fourier Transform as:

$$F(p)=\int_{-\infty}^{\infty}e^{-pt}f(t)dt$$
$$F(\omega)=\int_{-\infty}^{\infty}e^{-j\omega t}f(t)dt$$

and stability indicates that the real part of $p$ must lie between $a$ an $b$.(For stability I'm guessing?)

As per the textbook, if we take the values of $p$ to be purely real and ignore the imaginary axis, then $F(\omega)$ doesn't even exist.

logically, I can't see how that could happen.
I was wondering if anyone could explain how that is so.
As well as why we limit the real part of $p$.(guessing it's similar to the unit circle in Z transform).


texts: Extended Kalman Filter - how do I get transition functions

texts:
I am measuring position and velocity, both have some noise in them. Velocity is defined as derivative of position. The system is apparently non-linear so I need to use EKF.

Model:

\begin{align}
X&=\begin{bmatrix}p\\v\end{bmatrix}\\
p_{k+1}&=A_0\sin\left(b_0k\right)+A_1\sin\left(b_1k\right)+A_2\sin\left(b_2k\right)\\
v_{k+1}&=A_0b_0\cos\left(b_0k\right)+A_1b_1\cos\left(b_1k\right)+A_2b_2\cos\left(b_2k\right)\\
\end{align}

Questions:

I suppose this is a non-linear system. Am I right?
If yes, I should use Extended Kalman filter. But how do I get $f$ and $h$ functions for my case (below)? The current state doesn't depend on previous state; except the $k$ (step) variable.
\begin{align}
x_k&=f\left(x_{k-1}\right)+w_k\\
z_k&=h\left(x_k\right)+v_k
\end{align}


texts: Restrict a state value in an Extended Kalman Filter

texts:
I am trying to use an EKF to filter/predict a noisy sinusoidal waveform using the states of phase, rate of change of phase and amplitude. I am finding, however, that sometimes the EKF amplitude state estimate is a negative value which, in the context of my problem does not make any sense, although a negative amplitude is mathematically plausible as it can represent a phase shift that flips or inverts the waveform.

Is there any way to limit the state evaluations within the EKF such that the amplitude will always be a positive value? This will, of course, also have implications for the estimation of the phase state?


texts:
Edit: Could the following be the answer?

Generate WGN-like-signal which is centered around a set dBm value.
Treat that signal like it was a frequency domain representation of an unknown X time-domain signal.
Do some math to get he unit back into Volts/Hz from dBm/Hz (using reference impedance).
Run inverse FFT.

Original Question:

I am looking for some advise on how to generate an array of samples pulled from normal (or Gaussian) distribution.

But, and here's the deal.
I want to control its PSD or ASD "directly", assuming I am using a refernce impedence of 50 Ohms.

I would like to know Ahead-Of-Time (AOT), before signal gets generated, what STD Dev should be applied to achieve a set PSD or at least average spectral level in a given bandwidth/channel.

This value should be set by a user.

Code below should prove I know (more or less) what PSD or ASD is, I am just looking for an algorithm advise, maybe someone did something similar?

Here's how I calculate the PSD or ASD for a generated white-noise-like signal in a standard way using numpy.

import numpy as np
import matplotlib.pyplot as plt
import logging

if __name__ == "__main__":
    logger = logging.getLogger("Main")
    logging.basicConfig(format='%(levelname)s:%(message)s', level=logging.DEBUG)
    DRAW_PLOTS = False
    mu, sigma = 0, 0.001  # mean and std deviation
    N = 1000
    timestamp = 1  # 10 Hz sampling rate

    # This is in Volts and no timestamp
    s = np.random.normal(mu, sigma, N)
    # Because signal is in Volts
    # our result of fft is in Volts/Hz or would be if we say our timestamp for
    # signal `s` is eqault to 1 second
    # This noise_spectrum can be called Amplitude Spectral Density (ASD)
    # If we want PSD we have to go from ASD to PSD:
    # Square the voltage and also maybe apply some reference impedance
    # ... or don't if you want to end up with Voltz (^2) / Hz
    noise_spectrum = np.fft.fft(s)
    # The following is in dBV^2/Hz
    noise_spectrum_db = np.multiply(
        20,
        np.log10(noise_spectrum)
    )
    # The following is in dBW/Hz (divided by 50 Ohm)
    noise_spectrum_dbw = np.multiply(
        20,
        np.log10(noise_spectrum/50)
    )
    freq = np.fft.fftfreq(s.size, timestamp)
    logger.info(
        f"Avg channel power in dBW: {np.real(np.mean(noise_spectrum_dbw))}"
    )
    logger.info(
        f"Avg channel power in dBm: {np.real(np.mean(noise_spectrum_dbm))}"
    )
```


texts: Help me understand how Pulse-Amplitude Modulation (PAM) works

texts:
Reading on the Internet i learned that PAM is basically a sampling technique which converts an analog signal into a discrete signal. It modulates a train of pulses by varying the amplitudes of the individual pulses in a regularly timed sequence based on the analog wave. (Multiplying the analog signal with a train of pulses)

In my mind i have the idea of PAM being similar to AM but instead of an analog carrier we have a train of pulses.

On the other hand in some textbook i read that the way PAM works is by grouping bits from a stream into symbols and each symbol corresponds to a pulse with a certain amplitude.

This seems to me as a completely different process than the first definition.
So what is true? and if both are true when is the second way to do PAM used?


texts: Power Spectral Density as a single number confusion

texts:
I'm trying to recreate the results of a machine learning applied to the DSP classification in the article: link.
I have a signal (activity measurements from a smartwatch) per patient, so about 30 signals in total, quite long. For machine learning algorithm authors take traditional approach: extract statistical values about the signal, create columns (1 column per feature, 1 row per patient/signal) and plug this into a classifier. I'm confused about the spectral density used as feature.
Questions:

If I take PSD of a signal, I get another signal (time series) with values, not a single value, am I right? Therefore I would get not 1 column with spectral density feature per patient, but a full signal.

Suppose the article is imprecise (it is in other places) and authors made some "mental shortcut" while writing this. What else could they mean by "spectral density" feature? It should be 1 column, so 1 number extracted from each signal from patient. They mention it is in frequency domain, so PSD would be right, but it would not produce a single number.


texts:
I am reading "The Scientist and Engineer's Guide to Digital Signal Processing" and trying to understand Figure 3.5 below which is about the sampling theorem and aliasing.

I do not understand the picture of the frequency domain. The signal can be represented as
$$f(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}F(x)e^{ixt}dx$$
where $F(x)$ is the Fourier transform of $f$. This represents $f$ with the "frequencies" $e^{ixt}$ and these can have $x$ negative. How then am I suppose to understand Figure 3-5 which shows the frequency domain for only greater than zero and moreover real and not complex.

Also why does the frequency domain of the impulse train have the original frequency spectrum appear flipped half the time? My current understanding is that the repetition comes from something like
$$\mathcal{F^{-1}} \{f(t)\Delta(t)\} = F(x)*\Delta(x)$$
where $\Delta$ is the dirac comb and probably some kind of constant like $2\pi$ is needed in there somewhere. How does the flipping half the time in the frequency domain arise from $F(x)*\Delta(x)$. The image I think for convolution is that $\Delta$ is flipped and as it is moved to the right by $s$ the "area" of $F(x)\Delta(s-x)$ is the value of the convolution at $s$. This would seem to replicate the original spectrum multiple times but without the flipping half the time as shown in Figure 3.5.


texts:
It seems to me intuitively if I have a signal $y$ with sampling frequency $f_0$, then if I resample that signal with $f_{new}$ the correct number of samples will be $\frac{f_{new}}{f_0} \cdot len(y)$.

But I'm using a library which includes an example (line 26) where the length of the new signal is calculated as $\frac{f_{new}}{f_0 +f_{new}} \cdot len(y)$. It seems to me this is a mistake? If I have a signal with 48000 samples at the same $f_s$ and resample it with $f=16000Hz$ then the new signal should have length 16000 right?


texts:
I'm trying to find out if the correction (Jackson, Nelatury, Mecklenbräuker) could improve the (IIM based) filter response near Nyqvist.

Here's my c++ routine (which is used for to calculate various RIAA and non-RIAA filters by just changing the timeconstant values and samplerate):

double a0, a1, a2, b0, b1, b2;
double fs = 44100;
//timeconstants (case RIAA):
// frequency -> time conversion 1/(2*pi*fc) (= R*C)
//poles
double p1 = 3180e-6; // 1/(2*pi*50.05Hz)
double p2 = 75e-6;    // 2212Hz
//zeros
double z1 = 318e-6;  // 500.5Hz
double z2 = 0.0;       // 3.18e-6 for Neumann pole (50kHz)

double pole1= exp(-1.0/(fs*p1));
double pole2 = exp(-1.0/(fs*p2));
double zero1 = exp(-1.0/(fs*z1));
double zero2 = exp(-1.0/(fs*z2));

a0 = 1.0;   // = 1.0
a1 = -pole1 - pole2; // = -0.931176
a2 = pole1 * pole2; // = 0
b0 = 1.0; // = 1.0
b1 = -zero1 - zero2; // = -1.731986
b2 = zero1 * zero2; // = 0.733838

Tried to google "ready to use" solution of this but the only source code I found few papers I could use for correction.

Papers:

Jackson
Mecklenbräuker
Nelatury

Juha

EDIT:

Bypassed the above mentioned papers. Improved the method so that now there's no need for additional correction biquad but by using z2 (which is unused) for the correction. One can decide where the error lies. Result:


texts: STFT on time varying signal with good time and frequency resolution

texts:
I am trying to determine the main frequency of a noisy signal that varies in frequency over time. Ideally I want to detect changes in the frequency as rapidly as possible - say 50Hz update rate, but I also want to get as close to the actual signal frequency as possible. The signal I am looking for varies between 100Hz and 300Hz and I am sampling the signal at 1Khz.

Currently I am using a STFT of length N applied every N/2 samples (so hop size N/2) and using a Hanning window on the input. I am then averaging the 3 overlapping frames via Welch's method to achieve a final estimate of the signal I am looking for.

I am aware that the maximum frequency resolution I can get is 1000/N so I am using Jain's method to interpolate between bins. I am also aware that the maximum time resolution I can get is 1000/N Hz. I know that these two are in tension.

Sometimes the signal is missing. In order to detect this I am calculating the standard deviation as well as the average of the three overlapping frames and using the ratio to determine the "noisiness"

My questions are:

Is this a valid procedure?
Is using Welch in this context valid given that the input signal varies?
Am I averaging the right data?
Is the time resolution really 1000/N or is the fact that I can produce output at 1000/2N Hz better?
Would more overlap help? If so how?
Is my noise calculation valid?

I have heavily relied on https://holometer.fnal.gov/GH_FFT.pdf


texts: Accuracy of IFFT(FFT(X)) unit transformation

texts: I want to use FFTs to do seriously accurate interpolation on band-limited data.  To do that, I need to get a handle on the fundamental accuracy of the FFT() and IFFT() algorithms available.  My idea is to look at IFFT(FFT(X)) which should be a unit transformation.  I want to get the maximum difference between the input and output data down below the region of 1 part in 10^8 or even 10^9.  Currently I'm achieving ~1 part in 10^6 using 64-bit float, but I don't know if that is in line with what is to be expected.  I'm thinking that surely someone has done all this before.  Does anybody know the answers?  What can expect if I try using higher-precision data types?  With 64-bit float speed is not a concern.


texts:
I am undergoing a course on DSP and we were recently introduced to the amazing concept of DFT.

But I have some apprehensions against it. DFT is supposed to be the methodology of fourier transform, as implemented in hardware or software. Now, if I want to convolve two signals' DFT, I could simply linearly convolve them. There would be no wraparound error, because all that I have to do to implement linear convolution is tell my computer that the array now flips and then multiplies and slides - i.e., linearly convolves.

Circular convolution was developed to avoid wraparound error; however, the computer does not know that DTFT is a periodic function, and, thus, shifting the DFT array will render me with a circular shift.

I know my opinion is rather radical, but I am sure that this is because I have a very abstract level acquaintance to this concept. Can someone enlighten me more, herein?


texts: Is it possible to implement FIR band pass filter using two FIR band stop filters?

texts: I am trying to filter a signal using a given FIR band stop filter. I have to implement a band pass filter.  Assuming I want to get the signal between normalized cutoff frequencies [0.2,0.3] - can I use two band stop filters over the same signal - one for the range [0,0.2) and second for the range (0.3,1] ?


texts:
I am designing an IIR Filter with the following designs constraints:

Passband $(0.75\pi,\pi)$
Stopband $(0,0.70\pi)$
No more than $0.5$ % ripple in the passband
No less than $40$ dB attenuation in the stopband

Using the designfilt command I am able to design the filter.

However when I am using the following code for the "Butteworth" filter I do not get any Phase response. The magnitude response is fine, but the Phase values have NaN values - "Not a Number". The same codes are able to design for "Chebyshev" and "Elliptic" filters.

Wp = 0.75;
Ws = 0.7;
Rp = 20*log10(1.005);
Rs = 40;

[Nb, Wnb] = buttord(Wp, Ws, Rp, Rs); % Nb: minimum order; Wnb: cutoff frequency
[Zb, Pb, Kb] = butter(Nb, Wnb,'high'); % Zb : zeros ; Pb : Poles , Kb : gain

SOSb = zp2sos(Zb, Pb, Kb);

freqz(SOSb)

Response for Butterworth Filter

Response for Chebyshev I filter


texts: Constellation diagram for QPSK modulation

texts:
I need to plot (in MATLAB R2011a) constellation diagram of a QPSK modulated signal after AWGN was added to the signal.  My current MATLAB code looks like this:

b=round(rand(1, Nbits));
fill=10;
T=1/1000;
fs=120000;
fc=4000;
t=0:2/fs:length(b)*T-1/fs;

for k=1:floor(Nbits/2) %Dividing pair and odd bits
    even(k)=b(2*k);
end
for k=1:round(Nbits/2)
    odd(k)=b(2*k-1);
end
E=A*sqrt(T/2);
for k=1:length(even)  %NRZ code
    if (even(k) == 1)
        evenNRZ=[evenNRZ ones(1, fill)*E];
    else
        evenNRZ=[evenNRZ -ones(1, fill)*E];
    end
end

for k=1:length(odd)
    if (odd(k) == 1)
        oddNRZ=[oddNRZ ones(1, fill)*E];
    else
        oddNRZ=[oddNRZ -ones(1, fill)*E];
    end
end
fi1=sqrt(2/Ts).*cos(2*pi*fc*t); %Basic functions
fi2=sqrt(2/Ts).*sin(2*pi*fc*t);

up=fi1.*evenNRZ;
down=fi2.*oddNRZ;

QPSK=up+down; %QPSK signal

QPSK_n=awgn(QPSK,20,'measured'); %QPSK signal + noise

This is the method we are required to use to make QPSK modulation simulation for an assignment. Now I need to plot constellation diagram from QPSK_n signal. Anyone has any idea on how to do it?


texts:
This can be a very silly question, but I'm quite confused:

If we take the Fourier transform of any constant signal, we get an impulse at zero, which says that its frequency is zero and, hence, it is non-repeating and its period is infinity.
By the definition of a periodic signal, $$F(z+p) = F(z)$$ therefore, for a constant function $$F(z) = c$$ $$F(z+p) = c, \forall p >0$$ Therefore, a constant function is periodic, but its period is undefined or can be defined as anything.

Which of the above arguments are correct and why? Please clarify. Thanks in advance.


texts: Can we see the Fourier transform as a filtering operation?

texts:
Is the Fourier Transform (FT) a filtering operation?
Is it possible to graphically represent this transform?
I know that a signal can be represented in the frequency domain, but I want to know if the FT operation can be "drawn" as a coefficient series.

The FT matrix is composed by $N\times N$ elements, and I can apply it to the input sequence of length $N$ in order to obtain an output sequence of length $N$. But I want to understand if there is another way to see the FT matrix in one dimension.


texts: How is signal to noise ratio actually measured by receiver equipment?

texts:
This sounds like quite a basic question but it surprised me, how is SNR actually measured?
You have the incoming signal:

It seems like the SNR is just the visual comparison of the peak signal 'amplitude' to the noise 'amplitude'. In this case the signal 'amplitude' clearly includes noise as well, because it is summed together i.e. it is the sum of the noise 'amplitude' at that frequency and the signal 'amplitude' at that frequency because 'amplitude' in this case referring to volt-seconds under the curve of a complex sinusoid at that frequency times the signal (which contains the noise sinusoid and signal sinusoid at that frequency, which creates a greater amplitude sinusoid at that frequency and therefore there is more area (volt-seconds) under the curve when that sinusoid is multiplied by a complex sinusoid of the same frequency.
I'm assuming the received signal power would be calculated by squaring the $V_{RMS}$ of the samples taken and I assume noise is also sampled at a separate occasion to get the $N_{RMS}$ (or so called noise variance, but real noise doesn't have a 0 mean i.e. no 0 frequency component), or it uses a frequency spectrum analyser or something to get the voltseconds at unused frequencies and then integrate it with respect to frequency i.e. $s^{-1}$ to get the average voltage and then square it to get the average power, the problem with that being that you can't measure the noise on the frequencies that are used, so I would think that it is acquired separately, like in the wifi layer 1 preamble and training sequences.
But the signal comes in as the signal+noise, so surely it's actually (signal (${V_{RMS}}^2$ + noise ${V_{RMS}}^2$) / noise ${V_{RMS}}^2$) that is being measured, unless noise is deliberately subtracted from the signal i.e. received signal (${V_{RMS}}^2$ - noise ${V_{RMS}}^2$) / noise ${V_{RMS}}^2$ to get the SNR? I'm assuming that the traditional $P_{\text{signal}}$ means the signal without the noise and not the received signal as it is seen.


texts: Determine filter given transfer function

texts:
I am given $$H(z) = 1 + \frac{\alpha}{1-\alpha z^{-1}}$$ where alpha is between $0$ and $1$. This is apparently is a low-pass filter with cutoff frequency $f_c$.

How can I see that? And how can I compute/plot the frequency response in Matlab? Neither freqz(1+alfa,[1 -alfa]) nor abs(H(z)) seems to yield a low-pass filter on the above transfer function in Matlab.

Furthermore, $1-H(z)$ should be a high-pass filter. But again, using freqz(-alfa,[1 -alfa]) in matlab it is clearly not a high-pass filter.


texts:
I everyone,

I'm a computer scientist but not from the signal processing community. I want to implement a real time digital low pass filter (to begin with) for a sound processing application. I've done some research and found some solutions:

Wikipedia gave me this filter: y[i] := y[i-1] + α * (x[i] - y[i-1]),
which is exactly the kind of formula I want.
I've found an online filter generator, that can generate C code from the description of a Butterworth filter. Great.

So my question is not really "give me a filter algorithm". Its more:
How do you make those formula?

What I want is not a ready-to-use algorithm that I can copypast. I want to understand what I'm doing and how it works. However! All the documentation that I've found until now, including the Wikipedia pages, are directed to people that already have a solid background in signal processing (a lot of unexplained notations, function names, etc.). This is not my case.

I need someone to explain to me how a low pass filter works, how the Butterworth filter (or something else I don't care) works, how I can extract a nice algorithm from it with the correct coefficients. And all that like I'm a 5 years old.

I haven't found anything like that yet. If you'r willing do it here you'r welcome. Or you can give me links to online literature, blog posts, anything that popularize all thoses notions. You'll make me happy.

Thanks <3

EDIT: I want to clarify that I'm not looking to know how a general low pass filter works, using a convolution and all that. That I already know. I mentioned the Butterworth filter and associated algorithm because it is a linear filter (algorithmically speaking). Something that I can reasonably use in a real-life application for real time processing.


texts: Principal and Actual Values for Phase Response

texts:
Is there  a way to retrieve both Principal and Actual values for a Phase Response using the fvtool (Open Filter Visualization Tool) in MatLab?

Thanks


texts: Scale and Rotation invariant color based feature descriptors

texts:
SIFT/SURF/ORB etc are grayscale based, if the two grayscale bitmap is very similar, but the RGB channel is quite different, how to classify them? So, I think there should be some color based feature descriptors, also have scale and rotation invariant, like SIFT/SURF/ORB, but for RGB or HSV.

But util now, I cannot found such a feature descriptor ......


texts: which spatio-temporal properties does dynamic textures exhibits?

texts:
know that dynamic textures are spatially repetitive and temporally varying patterns.

1.I read one sentence in one paper i.e. They are sequences of frames exhibiting certain stationarity properties over time. They exhibit spatial self-similarity and temporal self-similarity. but i don't understand how it does?

2.Also how examples like "traffic on road" or "rafting" can be considered as Dynamic Textures? Because i don't thing they exhibit any stationarity property over time.is it?


texts: Is it meaningful to find linear convolution of just two random sequences?

texts: Could it bring any meaning? Is this kind of convolution be useful is solving anything?