texts: Removing periodic noise from image using Fourier transform

texts:
I'm trying to get rid of some periodic flicker noise through post-processing of the recorded images. The reason for these artifacts is that the electronic rolling shutter of the camera reads each line sequentially, and this mechanism results in brightness variations in the image, when the picture is taken under fluorescent lighting (due to 50 or 60 Hz AC power frequency).

I tried calculating the Fourier transform of the image and suppressing the components that might cause flickering. I inspected the magnitude spectrum (as you can see below for an example image) and set the vertical components around the DC component to a very small value (I preserved the DC component). In the end, I managed to eliminate the flicker because I suppressed the flicker frequency, but also the image was sort of deformed since I, probably, also suppressed some other frequencies that contributed to the signal's energy.

I would like to know what would be a better way to detect the flicker frequency and eliminate/suppress only that frequency? In my current implementation, I could be doing something fundamentally wrong, so I would appreciate any guidance and hints. Below, I show the resulting image that I got for the kind of filtering that I explained above. You can also find my Python code below.

import cv2
import numpy as np
from matplotlib import pyplot as plt

img = cv2.imread('images/flicker2.jpg',0)
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
# calculate amplitude spectrum
mag_spec = 20*np.log(np.abs(fshift))

r = f.shape[0]/2        # number of rows/2
c = f.shape[1]/2        # number of columns/2
p = 3
n = 1                   # to suppress all except for the DC component
fshift2 = np.copy(fshift)

# suppress upper part
fshift2[0:r-n , c-p:c+p] = 0.001
# suppress lower part
fshift2[r+n:r+r, c-p:c+p] = 0.001
# calculate new amplitude spectrum
mag_spec2 = 20*np.log(np.abs(fshift2))
inv_fshift = np.fft.ifftshift(fshift2)
# reconstruct image
img_recon = np.real(np.fft.ifft2(inv_fshift))

plt.subplot(131),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(132),plt.imshow(mag_spec, cmap = 'gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.subplot(133),plt.imshow(mag_spec2, cmap = 'gray')
plt.title('Magnitude Spectrum after suppression'), plt.xticks([]), plt.yticks([])
plt.show()

plt.figure()
plt.subplot(121),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(img_recon, cmap = 'gray')
plt.title('Output Image'), plt.xticks([]), plt.yticks([])
plt.show()


texts: 90 Degree Phase Shift / Rotation Algorithm for SDR

texts:
I'm working on a library to preform frequency demodulation from an SDR signal. I'm using https://github.com/sysrun/rtl-sdr/blob/master/src/rtl_fm.c as a crib for this. The algorithm starts by taking the samples (I & Q data in as 8-bit unsigned values I believe) from the radio and putting them through a 90-degree rotation. From this answer to a related question, I believe I understand why this is needed, but I don't understand how the algorithm achieves this. The algorithm is as follows, where buf holds the (I & Q) samples directly from the radio:

void rotate_90(unsigned char *buf, uint32_t len)
/* 90 rotation is 1+0j, 0+1j, -1+0j, 0-1j
   or [0, 1, -3, 2, -4, -5, 7, -6] */
{
    uint32_t i;
    unsigned char tmp;
    for (i=0; i<len; i+=8) {
        /* uint8_t negation = 255 - x */
        tmp = 255 - buf[i+3];
        buf[i+3] = buf[i+2];
        buf[i+2] = tmp;

        buf[i+4] = 255 - buf[i+4];
        buf[i+5] = 255 - buf[i+5];

        tmp = 255 - buf[i+6];
        buf[i+6] = buf[i+7];
        buf[i+7] = tmp;
    }
}

The comment is somewhat confusing, as I don't understand how the complex values match the array. I assume the values in buf correspond to the parts of the complex numbers like such: buf[i] + buf[i+1]j. If that's the case, then why is this a 90-degree rotation? Why does it require 4 complex values?

(buf[i+0] + buf[i+1]j) -> ( buf[i+0] + buf[i+1]j)
(buf[i+2] + buf[i+3]j) -> (-buf[i+3] + buf[i+2]j)
(buf[i+4] + buf[i+5]j) -> (-buf[i+4] - buf[i+5]j)
(buf[i+6] + buf[i+7]j) -> ( buf[i+7] - buf[i+6]j)

Any insight into this algorithm would be greatly appreciated... thanks!


texts:
Can anyone explain what spectral entropy is?

Does noise with a restricted bandwidth have the same spectral entropy as white noise?


texts: Filtering $ n \times n $ Images by Separable $ m \times m $ Filters. Computation Time for Filtering Using FFT, 2D Convolution and Two 1D Convolutions

texts:
Consider filtering square $n\times n$ images by square, separable $m\times m$ filters.

Give general equations for the computation time for the following approaches to this:

Filtering using FFT
2-D convolution
Two 1-D convolutions

Use capital letters to indicate implementation-dependent constants.

For particular implementations of the methods of filtering in part 1) the following computation timings were recorded:
n=64    m=3   FFT=2,498,561      2-D convolution=73,829      two 1-D convolutions=245,761
n=256   m=5   FFT=53,084,161     2-D convolution=1,179,649   two 1-D convolutions=3,932,161
n=1024  m=3   FFT=1,059,061,761  2-D convolution=18,874,369  two 1-D convolutions=62,914,561
n=1024  m=11  FFT=1,059,061,761  2-D convolution=471,859,201 two 1-D convolutions=314,572,801

Estimate the values of the implementation-dependent constants to 1)
For the implementation in 2) give values of $m$ and $n$ for which the FFT implementation is expected to be fastest.

My point of view about:

FFT is: Firstly, applying FFT to each column separately and then to each row. FFT for $N$ data points takes about $2N\log N$ multiplications, for $n\times n$ images it will be $2n\cdot 2n\cdot\log n$ multiplications. Secondly, inverse FFT to get back to sensible image. Another $2n\cdot 2n\cdot \log n$. Thirdly, multiply by FT of the kernel, it will be $n^2$ complex multiplication, that is $4\cdot n^2$ real multiplication. Thus the computation time for FFT is $\mathcal O_1\left(8\cdot n^2\log n +4\cdot n^2\right)$. Is that correct?
2-D convolution: I am confused about this. It could be $\mathcal O_2\left(n^2\cdot m^2\right)$. Is it correct?
Two 1-D convolution: Firstly, convolution per row. It should be $(n+m-1)^2$. Secondly, transpose the matrix of the image. Will this process include the computation time? Thirdly, convolution on the rows again. Another $(n+m-1)^2$. Fourthly, second transpose, return the image back to its original orientation. Computation equation is $\mathcal O_3\left(2\cdot\left(n+m-1\right)^2\right)$.

However, according to my computation equations, I cannot compromise my equations to the recorded computation timings. In this way, I think I need to correct my equations to solve the problems in part 3).


texts: What is the technical term to join two signals with different phases?

texts:
there is a name for this:
If I have two signals with a discontinuity, arrange the phases of one signal to fit the other and not see the discontinuity.


texts: How do I optimize the window lengths in STFT?

texts: I have many EEG signals and I want to analyze them using linear methods such as STFT (Short Time Fourier Transform). In STFT , How can I optimize the analysis window length, to reflect the frequency spectrum of each analysis window in a proper way?


texts:
i am studying the methodes CF minimization. and  have a doubt.

in each methods there are steps with DFT and IDFT or like this FFT and IFFT.
Everyone provides the fact, that they use DFT and IDFT for calculation muti-sin with new phase and then take DFT or FFT to present in time domain.

But what is the purpose of dft and idft in minimization of the CF?
Are their function  to represent the signal in F -or T- domain?


texts:
One way to describe a practical integrator ("leaky integrator") is

$$
H(s) = \frac{g R}{1 + sRC}
$$

I am trying to understand how to choose the values $g$, $R$ and $C$ because in practice, I will have constraints on them.

Here is my attempt: Suppose for a moment, that $R \rightarrow \infty$, then it becomes an ideal integrator

$$
H(s) = \frac{g}{sC}
$$

When I just look at this equation (or look it its time domain version) it seems that g and C do not matter because it's just a scalar. Do I understand it correctly that $g/C$ relates to saturation (practically, there is a maximum output value). Suppose I want to implement a finite integration:

$$
x_i = \frac{g}{C} \int_0^{T_w} x(t) \, dt
$$

and I know that $x(t)$ and $x_i$ can have a maximum value of $K$.
Then I could derive (to avoid saturation in the worst case):

$$
K < \frac{g}{C} T_w K \Rightarrow \frac{g}{C} > T_w
$$

This would restrict the ratio between $g$ and $C$. Does this make sense?

Back to the non-ideal integrator used for a finite integration interval. In this case, $1/RC$ can be interpreted as the first pole and after this frequency, the ideal integrator starts.

Since I integrate over a finite period $T_w$, is it valid to say that the lowest integrateable frequency is $1/T_w$ and hence $1/RC \ll 2\pi/T_w$ ?


texts: Laplacian of Gaussian Approximation and Gaussian Blur as the Solution of Heat Equation

texts:
While I was reading SIFT paper(Lowe, 2004), I came across the method that he apply "heat diffusion equation" to Gaussian function to derive that
$$ \frac{∂G}{∂σ}
= σ∇^2G
$$
I searched Wikipedia and found that original heat equation is

$$ \frac{∂u}{∂t}
= σ∇^2u
$$
I wonder how can time derivative equation be applied to sigma derivative equation. Is there anybody to understand me?


texts: Why there is Difference between shapes of ROC of z domain and s domain?

texts:
ROC(region of convergence) of Z domain is shown by a circular region while ROC in S domain is shown by a rectangular(approximately looking like rectangle) region

What is the reason of this difference in shapes of ROC regions?


texts: Formulas to plot Frequency Impulse Response Coefficients

texts:
I made this question https://dsp.stackexchange.com/questions/19731/newbieplotting-fft-and-impulse-response-coefficients-using-java

I have the IRC: Impulse Response Coefficients:

$h(k), k = 0, 1,.., N-1$

Low Pass: $h_D(k)=\sin(k\omega_c)/{k\pi} , k\neq0; h_D(0)=2f_c$

High Pass: $h_D(k)=-\sin(k\omega_c)/{k\pi}, k\neq0; h_D(0)=1-2f_c$

Band Pass: $h_D(k)=(\sin(k\omega_{cH})-\sin(k\omega_{cL}))/{k\pi}, k\neq0; h_D(0)=2(f_{cH}-f_{cL}) $

Band Stop: $h_D(k)=(\sin(k\omega_{cL})-\sin(k\omega_{cH}))/{k\pi}, k\neq0; h_D(0)=1-2(f_{cH}-f_{cL}) $

Hamming Window Function:

$w(k)=0.54+0.46\cos( \frac{2\pi k}{N} )$

so we have:
$h(k) = w(k)\times h_D(k)$

obviusly my $h(k)$ can be represented as:

double[] h = new double[N];  //N:odd integer
...
//hD[(N-1)/2] -> center or hD(0)
...
for (int  k = 0; k < N; k ++) {
  h[k] = hD[k]*w[k];
}

I want to get the formula to evaluate the h[k]

In matlab to do this:

[hFlt,fFlt] = freqz(h,1,256);
mLP = 20*log10(abs(hFlt));
plot(fFlt,mLP);

I want to obtain hLP and fLP values to plot filters response (without using matlab).

But I don't know the formulas(procedure or algorythm) to plot the frequency response.


texts: Inverse Continuous Wavelet Transform derivation?

texts:
Wiki writes iCWT as
$$
f(t) = C_{\psi}^{-1} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} W_f(a,b) \frac{1}{|a|^{1/2}} \tilde\psi \left(\frac{t - b}{a}\right) db \frac{da}{a^2}, \tag{1}
$$
where $\tilde\psi$ is the dual wavelet of $\psi$ in CWT, $W_f$ is the CWT of $f$, and $C_\psi$ is the admissible constant
$$
C_\psi = \int_{-\infty}^{\infty} \frac{\bar{\hat{\psi}}(\omega) \hat{\tilde{\psi}}(\omega)}{|\omega|} d\omega, \tag{2}
$$
Hat = $\mathcal{F}$, overbar = complex conjugate. Some wavelets are own duals, $\psi = \tilde\psi$; then the integrand is $|\hat\psi(\omega)|^2 / |\omega|$.

How does $(1)$ invert the CWT? And why is there a $a^{-2}$ normalizing factor? Above wavelets are L2-normed - do we normalize iCWT differently if they are L1-normed (i.e. $|a|^{-1}$)?
The normalization is particularly troubling; in forward-transform, we do $|a|^{-1/2}$, so in inversion one might expect $|a|^{1/2}$ for cancellation, but it only goes the other way: another $|a|^{-1/2}$, along $a^{-2}$. The role of $C_\psi ^{-1}$ is also unclear.


texts: Rotating and Shifting image doesn't change FFT2

texts:
I have an image and its fourier transform.
When I rotate it, its fourier transform rotates too, but I can't figure it out. Why does this happen?
On the other hand, when I shift the image, its fourier transform doesn't change. As I know, time shifting means frequency shifting. Am Iwrong?


texts: An alternative algorithm for blob detection

texts:
Here's a sample image of my experiment which is to be processed for finding the information of all the blobs (centre, dimension etc.). Basically the image contains some noise, few bold lines, and of course few blobs.

Full size image (12Mb)

My aim is to program a FPGA using an image processing package "Sapera APF". Sapera APF is an application for programming Teledyne DALSA frame grabbers to perform image processing. And unfortunately it doesn't allow us to write code in conventional way, neither allows to use basic coding features (if else statements, loops etc.) but only supports some inbuilt functions in form of code blocks, and we can only arrange these blocks.
FYI, Here's a link to a list of functions supported by Sapera:
https://www.dropbox.com/s/9f5qc9uki0lsalu/sapera_functions.docx

Now my first objective was to test my algorithm on matlab thereafter implementing it on Sapera. Here is my approach in matlab:

Applied mean filter (for step 3)
Thresholded image using adaptive thresholding (returning a binary image)
Removed lines by checking connectivity, all lines are connected with above 100 pixels      (for better connectivity mean filter was applied first)
Now since only blobs were remaining, their locations were identified (with grayvalue=255) and corresponding pixels were taken in original image.

This is a simple algorithm which works fine in matlab. But the real problem is that for step 2 (i.e. connectivity) of course I had to write loops going to each pixel checking for connectivity, which can't be done in Sapera. So I am looking for some other simpler algorithm to separate blobs in my images, without using loops and conditional statement, so that it can be implemented in Sapera, but now I just need its matlab code or algorithm.

Edit:
I cannot crop out strip containing blobs because other images may contain lines very close (maybe overlapping) to blob. For example, https://www.dropbox.com/s/m2j5pqo6stg4s1z/h125_s_-30.bmp


texts: How to perform 2D image rotation by FFT based image translation using matlab?

texts:
I am studying rotation using 3 translations as described in this scientific paper. According to the authors , it is possible to decompose a rotation into 3 shear translations which will be implemented in the frequency domain by FFT according to time-frequency shifting theorems.

For example the first translation is obtained as in equation (13) by doing the following operation for every pixel considering that the horizental coordinate x is variable and the vertical coordinate y is fixed:

Gx(x,y)=FFT-1(exp(-2*piau*y)*FFT(f(x,y)) where f is the value of the pixel x,y , a is a constant for the translation and FFT-1 is the inverse fourrier transform and u is x mapped into the frequency domain.

To code this transformation, one should think as follow:
We have to compute  Gx by changing the range of U from -M/(2*M) to (M/2*M) and by iterating y from 1 to N.

Surprisingly, the suggested implementation doesn't follow the procedure above:
Can somebody develop further why the implementation is as follow.

close all;
clear all;
I = double(imread('Put image here'))/255;
%imshow(I)
I=rgb2gray(I);
facteur=4;
[M, N, s] = size (I);

temp=ones(2*M,2*N);
temp(floor(M/2):floor(M/2)+M-1,floor(N/2):floor(N/2)+N-1)=I(:,:);
I=temp; clear temp;
imshow(I)
[M N s]=size(I)
 teta=pi/8;
a=tan(teta/2)
b=-sin(teta)

 Nx = ifftshift(-fix(M/2):ceil(M/2)-1);
 Ny = ifftshift(-fix(N/2):ceil(N/2)-1);
 Ix=zeros(length(Nx),length(Ny));
 for k=1:M
   Ix(k,:)=(ifft(fft(I(k,:)).*exp(-2*i*pi*(k-floor(M/2))*Ny*a/N)));
 end
 %imshow(Ix);


texts: Licence plate detection with different backgrounds

texts:
I just studied the Algorithm of ANPR (Automatic Number Plate Recognition) from the chp 5 of book Packtpub.Mastering.OpenCV.with.Practical.Computer.Vision.Projects in which they detect the vehicle licence plate , but their are some restriction in that algorithm , i call it restriction because its specific for one country and can not be used in my country , well in that algorithm they use the following algorithms to detect the licence plate area

Sobel filter
Threshold operation (otsu)
Close Morphological operation
Mask of one filled area
Svm classifier

But the thing they used is the background color of licence plate as their keypoint for the recognition and detection , but in my area licence plate have different backgrounds , most common are black , white and green , green is not much common but rarely used , black background used white number and white contain black numbers , they use floodfill algorithm for white background , but what should i use for 2 different backgrounds ? Size of the licence plate and size of character matter in any case ?

Here are the sample images of my regional vehicle's liecence plate


texts:
Polyphonic music transcription does not currently appear to be a solved problem.

How about the inverse of a small portion of the problem.  Are there any kind of spectral characteristics (from an STFT) that can be used to eliminate some musical chords from the probability space?  (e.g. this snippet of sound most likely does not contain any C# chord, or any kind of diminished minor chord, or this is a single note not a chord, etc.)

Assume the audio snippet is more-or-less stationary (transient attack removed, etc.), and that overtones for most or all individual notes are very likely present.  (And this question is not about inverted chords.)


texts: Finding n Amplitudes by DFT, what is correct normalization

texts:
Please forgive me if this has already been asked.
Let us assume an example with $x(t) = \sum_{i=1}^N A_i \sin(2\pi f_i t) $ under a given sampling frequency $f_s$, frequencies $\omega_i$ and Amplitudes $A_i$. My ideas was to calculate the DFT transform $\tilde{x}(\omega)$, find the peaks of $|\tilde{x}(\omega)|$ and infer the original amplitudes from them.

However, there are several obstacles.

There is leakage, and the Maximum of the peak is a bad indicator for the Amplitude of the signal; the normalization factor seems to be dependent on several factors.
The Peaks overlap.
I need a window function to get rid of the side lobes, which again changes my normalization.

Do you have experience in doing this? What is the best practice? Any useful literature?


texts:
I got a lab report to make in an image processing course in my university, and I was asked a question about the OpenCV camshift algorithm (for video tracking):

the camshift algorithm works on the hue component of the HSV color space.  I need to explain why working on this component should give better results than working in the RGB color space.

I must admit that I'm not familiar with OpenCV at all, and basically I'm leaning all the color model theory specifically for this course, so please consider that...thanks.


texts: Detection of number of yarns in the fabric image

texts:
I am a novice user in image processing and working on a project to determine number of yarns in a fabric image i.e. warps and wefts (Perpendicular and parallel yarns to x-axis).

I am trying to detect the yarns by taking the sum of columns and rows for warps and weft respectively but it doesn't seem working because the image is not very well focused and sharp. Also the method is much sensitive to the alignment of yarns in the fabric.

Please propose a better algorithm if possible for detection of yarns. Also if there is any possibility to make yarns straight for every captured image. I have tried IMROTATE but it creates the black regions in image which disturbs the sum of rows and algorithm gets failed.

Any help in this regard will be of much worth to me. Thanks.


texts: Answering Machine vs. Human -- Neural Network Features Selection

texts:
I've been tasked with creating an artificial neural network that can classify a telephone call as either answered by a human or an answering machine.
My knowledge of audio processing is, mildly put, limited. I'm quite new to machine learning as well, so feel free to question everything.  Here's what my run-of-the-mill backpropagation network looks like so far:
1. 37 input units;
2. One layer of 10 hidden units;
3. One output unit.

An audio file comes in, and I extract the data like so [X, fs] = wavread(filename). Then, the 37 features I obtain are the following:
1. The length of X;
2 - 21. In 5% increments, the mean of X in each interval (I know);
22. Spread;
23. The mean of short-time energy;
24. Zero crossing rate;
25 - 37. The mean of each in a set of 13 MFCCs.

Now, since I've obviously jumbled it all together, my question is the following: What do you think would be a good set of features for the problem at hand?

Thank you for your time


texts:
(From Shaums DSP outline, 2nd edition, page 248, problem 6.21)

Book says, evaluate the Sum:

$$
S = \sum^{N-1}_{n=0} \Bigg( x_1[n] \ x^{*}_2[n] \Bigg)
$$

when:

$$
\begin{aligned}
x_1[n] = \cos\left( \frac{2\pi n k_1}{N} \right) \\
\\
x_2[n] = \cos\left( \frac{2\pi n k_2}{N} \right)
\end{aligned}
$$

using Property:

$$
\sum^{N-1}_{n=0} \Bigg( x_1[n]\  x^{*}_2[n] \Bigg) = \frac{1}{N} \sum^{N-1}_{k=0}  \Bigg(  X_1[k]\ X_2^{*}[k]\Bigg)
$$

I start by converting cosines into complex exponentials:

$$
x_1[n] = \cos\left( \frac{2\pi n k_1}{N} \right)
$$

$$
x_1[n] =  0.5e^{j2\pi n k_1 / N}  + 0.5e^{-j 2\pi n k_1 / N}
$$

Applying Definition of DFT to find:
$$
X_1[k] = \sum_{n=0}^{N-1} \Bigg( 0.5e^{j2\pi n k_1 / N}  + 0.5e^{-j 2\pi n k_1 / N} \Bigg) e^{-j2\pi nk/N}
$$

$$
X_1[k] = \sum_{n=0}^{N-1} \Bigg( 0.5e^{-j(2\pi n/N)(k- k_1)}  + 0.5e^{-j(2\pi n/N)(k+ k_1)} \Bigg)
$$

From this I determine that sum over one period of a complex exponential is zero, except when k is selected to cancel out exponential function input to zero, such as when $k=k_1$ and $k=-k_1$.  Thus:

$$
X_1[k] = \begin{cases}
\frac{N}{2}&(k=k_1)\ or\ (k=-k_1) \\
0 & else
\end{cases}
$$

By similar means:
$$
X_2[k] = \begin{cases}
\frac{N}{2}&(k=k_2)\ or\ (k=-k_2) \\
0 & else
\end{cases}
$$

Here's where I get lost:

Books says, therefore:

$$
\sum^{N-1}_{n=0} x_{1}[n] x^{*}_{2}[n] = \frac{1}{N}\Bigg[ \frac{N^2}{4} + \frac{N^2}{4} \Bigg] = \frac{N}{2}
$$

Where does:
$$
 \frac{N^2}{4} + \frac{N^2}{4}
 $$
come from?


texts: Unexpected result after calculating the response of a system in frequency domain

texts:
I'm implementing a Python script that calculates the response of a system given an input. The system is $y(t) = x(t-2)$, i.e. it delays the signal by $2$, and the signal is $x(t) = \sin(3t)u(t)$.

As I'm learning the discrete Fourier transform (and it's inverse), I'm doing all the calculation in frequency domain. The procedure is very simple:

Sample $x(t)$ between $[-10, 10]$ where the interval between each sample is $0.01$
Calculate the DFT of $x(t)$ sampled (the result is $X(j\omega)$)
Multiply the $X(j\omega)$ by $e^{-j2\omega}$ to shift the signal (the result is $Y(j\omega)$)
Calculate the inverse DFT to obtain $y(t)$ sampled

In the end, the result is plotted and it's straightforward to check its correctness for $t > 0$. However, for $t < 0$, the output has some kind of strange behavior, as you can see below:

Am I missing something? Can someone please explain me what causes this?


texts: I Q sampling and baseband version of analytic signal

texts:
Is it correct to say that if we have a radio signal $s(t)$ centered around the angular frequency of $\omega$ as $\omega\pm\omega_B/2$ (where $\omega_B$ is the bandwidth of the signal) and the corresponding analytic signal is $s_a(t)=s(t)+j \hat{s}(t)$ (where $\hat{s}(t)$ is the Hilbert transform of $s(t)$) then $I$ and $Q$ as received from an SDR receiver (RTL-SDR) can be represented as

$I=Re(s_a(t) e^{-j \omega t})$
$Q=Im(s_a(t) e^{-j \omega t})$

In other words is it correct that $I$ and $Q$ are the real and imaginary parts of the analytic signal shifted to baseband?

I'm trying to derive the above from the assumption that the $I$ component is a low-pass filtered version of $s(t) \cos\omega t$ and the $Q$ component is a low-pass filtered version of $s(t) \sin\omega t$.

Considering that

$s(t) \cos\omega t = (s_a(t)-j\hat{s}(t))\times {1\over{2}}(e^{-j\omega t}+e^{j\omega t})$
$s(t) \sin\omega t = (s_a(t)-j\hat{s}(t))\times {1\over{2}}j(e^{-j\omega t}-e^{j\omega t})$

and then expanding these and discarding high frequency components (both negative and positive) I arrived at the following

$I={1\over{2}}(s_a(t) e^{-j\omega t} - j\hat{s}(t)e^{j\omega t})$
$Q={1\over{2}}j(s_a(t)e^{-j\omega t} + j\hat{s}(t)e^{j\omega t})$

And I expected that $I+jQ$ should give me the shifted analytic signal $s_a(t)e^{-j\omega t}$ but when I substituted $I$ and $Q$ from above into $I+jQ$ I ended up with $-j\hat{s}(t)e^{j\omega t}$ which is the negative frequency replica of the signal shifted to baseband, which would be a mirrored version of $s_a(t) e^{-j \omega t}$ (i.e. of what I expected).

Is my initial statement about the relation between $I$ and $Q$ and the analytic signal correct and there's just a mistake in my math or am I making a wrong assumption somewhere along the way?


texts: What is the advantage of applying a Hamming window for MFCC calculations?

texts:
I am trying to implement MFCC by following this and I found that there is a step in which signals are divided into frames and then passed on to Hamming windowing process. The reason given there is to correct the discontinuity at the start and last samples of the frame. OK that sounds good, then the output of the window is fed for FFT process. So I thought of doing this sample code separately.

My code:-

%Creating a signal to perform hamming windowing on.
fs=8000;
ts=1/fs;
N=1000;
t=(0:ts:(N-1)*ts);
x=sin(2*pi*60*t);

figure,subplot(221),plot(t,x,'.-'),title('Original Speech'),xlabel('Time');grid on
f=(-fs/2:fs/(N-1):fs/2);
subplot(223),plot(f,fftshift(abs(fft(x)))),title('Frequency Spectrum'),xlabel('Freq (Hz)');grid on
windowed=x.*hamming(length(x))';
subplot(222),plot(t,windowed,'.-'),title('Hamming Window Applied'),xlabel('Time');grid on
subplot(224),plot(f,fftshift(abs(fft(windowed)))),title('Frequency Spectrum After Hamming'),xlabel('Freq (Hz)');grid on

The plot generated is:-

My Problem is:-

What is the advantage of Hamming Windowing. I am a noob, so according to me, firstly we are loosing the amplitude at the staring and ending of the signal in the time domain and secondly, in the frequency domain, there is very little difference. So what is the advantage of doing all this much computation.
Please explain with as much detail as possible.

And it would be nice if someone could point out to more detailed (according to them) tutorial then what I am referencing or what they read while learning MFCC themselves.


texts:
PLL is typically used in carrier synchronization and its diagram is well known, consisting of a multiplier, a low pass filter, and a VCO.

In "Phase-locked Loop" at wiki (https://en.wikipedia.org/wiki/Phase-locked_loop), it says in the "Application" section that PLL can be used for symbol synchronization. Does anybody have a standard diagram for this application, with some analysis on how it works?


texts:
My question is whether the chirp z-transformed is capable of solving two very close tones that an FFT is not capable of solving. I understand that it is possible but I have done tests in Matlab and I cannot get more information from a signal.
The formulation that I have used is the following:$X_k= \sum_{n=0}^{N-1} x_n z_k^{-n}$ where $z_k=A W^{-k}$ and $W=e^{\frac{-j2\pi (f_2 - f_1)}{m F_s}}$ and $A=e^{\frac{j2\pi f_1}{F_s}}$ Fs is the frequency sampling, f2 and f1 are the start and final of the interest bandwidth and m is the number of samples.
Implementation
Fs = 8000;
L = 128;
T = 1/Fs;
t = (0:L-1)*T;
f1=480;
f2=520;
X = sin(2*pi*f1*t) + sin(2*pi*f2*t);
Y = fft(X);
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L;
figure
plot(f,P1)

f1=480;
f2=520;
X = sin(2*pi*f1*t) + sin(2*pi*f2*t);
m=length(X);
f1 = 100;
f2 =1000;
%Resolucion de 100-1000/128=7Hz
w = exp(-j*2*pi*(f2-f1)/(m*Fs));
a = exp(j*2*pi*f1/Fs);
Y2=czt(X,m,w,a);
fn = (0:m-1)'/m;
fz = (f2-f1)*fn + f1;

% Resolucion_czt2=fz(2)-fz(1);
hold on
plot(fz,2*abs(Y2/m))
axis([f1 f2 0 inf])
legend({'FFT', 'CZT'})


texts: How to reduce synthesis artifacts produced by phase vocoder?

texts: I implemented the phase vocoder algorithm in Python to time-stretch speech signals by following this paper and referring to this MATLAB tutorial. I can distinguish words in the original signal from the resynthesis result, and the pitch, as expected, doesn't get altered even though the time axis has been stretched(which is the goal of using phase vocoder). However, there exist noticeable artifacts. Specifically, the output contains more noise resembling Gaussian noise compared to the input, and the words somewhat sound like being spoken a robot. I have tried some improvements such as identity phase locking suggested by this paper, but the result remained unchanged. Has anyone encountered similar issues? How can I improve the result?


texts: What is issue in square wave reconstruction?

texts:
As shown in attached photo:

What is the reason that the approximate square wave $$x_N (t)$$ at right side is not like the original $$x(t)$$ at left side .?

Is the only reason is that number of coefficients N is very less(in this photo/case N=+_13 )and value of N needs to be raised ?


texts: Bandlimiting and OFDM spectrum - what is considered "bandlimited"?

texts:
I have been reading about OFDM modulation, it's not clear to me how you can have band-limited OFDM but given its prevalence in wireless communication there must be some designation of approximately band-limited that OFDM satisfies.

So my question is how is the bandwidth of OFDM systems calculated and if it's an approximation based upon the largest spectrum lobe, why is that good enough?

It doesn't even look like OFDM systems are pulse shaped to reduce the side-lobe strengths so what's going on?


texts: Determining the right order of time-domain low pass filter for desired cutoff frequency

texts:
I've been stuck on a problem regarding complex demodulation of a signal for a while and wondered if anyone could help.
Suppose I have a signal $X(t)$ such that
\begin{align}
    X(t) &= A(t)\cos(\omega t + P(t)) + \xi(t)  \\
    &= \frac{1}{2}A(t)[\exp(i(\omega t + P(t)) + \exp(-i(\omega t + P(t))] + \xi(t)
\end{align}
where the signal amplitude, $A(t)$, and phase, $P(t)$, are assumed to vary slowly relative to the frequency $\omega$ but be time-dependent.
The signal is demodulated via multiplication by $\exp(i\omega t)$
\begin{align}
    Y(t) &= X(t)\exp(i\omega t) \\
    &= \frac{1}{2}A(t)\exp(iP(t)) + \frac{1}{2}A(t)\exp(-i(2\omega t + P(t))) + \xi(t)\exp(-i\omega t)
\end{align}
We apply a low-pass filter to the demodulated signal to remove frequencies at or above $\omega$, thus removing the second and third term in the previous equation
\begin{equation}
\hat{Y}(t) = \frac{1}{2}\hat{A}(t)\exp(i\hat{P}(t))
\end{equation}
where $(\hat{\cdot})$ denotes a smoothed variable. We can finally extract the amplitude and phase with
\begin{equation}
    \hat{A}(t) = 2|\hat{Y}(t)|,\quad \exp(i\hat{P}(t)) = \frac{2\hat{Y}(t)}{\hat{A}(t)}
\end{equation}
For complex demodulation low-pass filtering can be achieved by convolving the signal with a window of length $N$. I wanted to test the quality of several window functions against a dummy signal, but to no avail. The main thing I needed to figure out is how to relate the filter order to the cutoff frequency in the low pass filter.
In my research I'm dealing with hour long signals and focusing on frequencies in the range 1-10mHz ($10^{-3}Hz$), so I make the dummy signal from the following:
def test_data(periods, noise=0, npts=1000, dt=1.0/24, phase = 0, amp=1):
    t = np.arange(npts, dtype=float) * dt

    x = noise * np.random.randn(npts)

    for p in periods:
        x += (np.cos(2 * np.pi * t / p + phase) * np.hanning(len(t)))

    return t, x
freqs = np.fft.rfftfreq(3600)*1000
freqs = freqs[(freqs>=1) & (freqs<=10)]/1000
t, x = test_data([1/f for f in [freqs[8],ulf[freqs]]], noise=0,
                     npts=3600,
                     dt=1)

Which yields the following

I believe that a good filter and/or filter order should be able to extract the correct amplitude information at the 2 frequencies that make up the signal, without much leakage into other frequencies. Since the cutoff frequency is at -3dB of the filter frequency response, I thought that seeking the order corresponding this attenuation would be the correct way forward:
from numpy.fft import fft, fftshift
res=[]
ff=[]

m = ['blackman','hanning','hamming','bartlett','kaiser','boxcar']
for method in m:
    rr=[]
    for nf in np.arange(51,3600,2):
        if method=='blackman':
            win = np.blackman(nf)
        elif method=='hanning':
            win = np.hanning(nf)
        elif method=='hamming':
            win = np.hamming(nf)
        elif method=='bartlett':
            win = np.bartlett(nf)
        elif method=='kaiser':
            win = np.kaiser(nf,3)
        elif method=='boxcar':
            win = boxcar(nf)

        A2 = fft(win,3601) / (nf/2)
        mag = 2*np.abs(fftshift(A2))
        freq = np.linspace(-0.5, 0.5, len(A2))
        with np.errstate(divide='ignore', invalid='ignore'):
            response = 20 * np.log10(mag)
        response = np.clip(response, -100, 100)
        rr.append(response)
    res.append(rr)

Yet running the complex demodulation at each $f$ with the corresponding -3dB order isn't performing as I'd expected

compared to setting $N=2/f$ which yields better results:

Am I just completely missing the mark here?


texts: In audio analysis, is the first rahmonic from the cepstrum equivalent to the first partial of a pitched sound?

texts:
I was reading this text from Oppenheim and Schafer, briefly explaining the concept of the cepstrum. I'm new to this concept.
He explains that the cepstrum is useful to quantize the distance between 'echoes' of a signal. This echo value can be obtained analysing the first peak in the “spectrum” of the log spectrum from the original signal, since it will bring with itself the time delay through the Spectral density and log operations.
He also describes this peak as a 'rahmonic'.
My question: considering an audio source like a voice or an musical instrument, is the first rahmonic related to the first partial of a sound (a.k.a. the pitch of a sound)?
I'm asking this because it seems these 'echoes', related to the timbre of a sound, are not only periodic, but the lowest partial in the signal.
Also, if both values are related, why a cepstrum representation, like MFCC, is useful? Aren't Frequency-domain representations (like STFT, CQT, chromagram, etc.) enough?


texts: expected value of two LTI output signals multiplied (cross correlation)

texts:
I have an input signal x (assumed to be iid Gaussian with $\mu=0$, $\sigma^2$) which is fed into two linear systems:

$y_1 = h_1 * x$
$y_2 = h_2 * x$

Now I would like to calculate $\mathbb{E}[y_1 y_2]$. What is the proper way to do this - preferably in the frequency domain?

Background: I want to calculate something like $\operatorname{var}(y_1 + y_2)$. This expands to $\operatorname{var}(y_1)+\operatorname{var}(y_2)+2\mathbb{E}[y_1 y_2]$. Now I am very familar with the first two components - the variance. I know that the variance of a stochastic signal is given by the autocorrelation sequence at position 0. With the Wiener-Kinchin theorem this translates to:

$$
\operatorname{var}(y_1) = r_{yy}(0) = \int_0^{\infty} \Phi_{yy}(f) df
$$

and then:

$$
\cdots = \sigma_x^2 \int_0^{\infty} |H_1(f)|^2 df
$$

Now for my problem - with $\mathbb{E}[y_1 y_2]$ I arrive at $\mathbb{E}[x^2 \cdots]$ (giving a variance) but the filter is not square magnitude.


texts: How to evaluate a background subtraction algorithm given ground truths?

texts:
I'm testing a background subtraction algorithm with sample image sequences. I am wondering how one can evaluate the accuracy of a background subtraction algorithm results given ground truths?

The only idea I have in mind right now is to take the difference between a result image and its corresponding ground truth image and calculate difference/total_white_pixels_of(ground_truth).

Cheers!

EDIT:
Dear all, by ground truth image I mean binary image with foreground as white blobs and background as black pixels. This is an example of input image and its ground truth:

source: http://wordpress-jodoin.dmi.usherb.ca/dataset2014/


texts: I am trying to do Noise cancellation using RLS algorithm at a software level, as part of my M.Tech Project work. Here, the error strength is continuously increasing, instead of decreasing as time grows. I need to get convergence after some time. Only thing i observed is that the correlation matrix values are increasing over time. Do I need to change parameters like forgetting factor, initialization of the correlation matrix or any other thing?


texts: Timing synchronization in limited channel

texts:
I need measure timing of an signal that is send through a noise band limited channel (wireless).
By measure timing, I mean measure as exactly as possible when the signal has arrived at the receiver.
My idea is to use a matched filter (correlator) to match some symbol send by the transmitter and then locate the maximum point to find the exact time of arrival.

But is this the best way to do it?


texts: Are there any template matching / optical flow algorithms for large image displacement?

texts:
I had some success with feature-based methods for image registration, but these methods are sometimes unsuitable (e.g. images containing repetitive patterns, very low contrast, low texture, no edges/corners).

I am therefore looking for a direct (pixel-based) global alignment alternative.

So far I looked on algorithms like Horn-Schunck, Lucas-Kanade, ECC but all the examples show only images with very small displacement (few pixels at most).

Is there any direct (pixel-based) method for registering images with a displacement of 200 pixels or more (in any direction)?

So far I am using simple template matching with exhaustive search over possible image shifts and Normalised Cross-Correlation error metric, which works, but it is extremely slow (basically $O(n^4)$) and does not allow for rotation.

Are there any faster methods allowing large image shifts/rotations?


texts:
I am trying to figure out a way to find out the vein pattern in the hand.. Please could someone suggest ways to do the same.? I know first i would need to find the ROI in the image and perform further processing .. Can gabor filter be used as one of the steps ?
Thanks


texts: Time delay estimation of oscilloscope signals using cross correlation

texts:
I have recorded 2 signals from an oscope. They look like this:

I want to measure the time delay between them in Matlab.
Each signal has 2000 samples with a sampling frequency of 2001000.5.

The data is in a csv file. This is what I have so far.
I erased the time data out of the csv file so that only the voltage levels are in the csv file.

x1 = csvread('C://scope1.csv');
x2 = csvread('C://scope2.csv');
cc = xcorr(x1,x2);
plot(cc);

This gives this result:

From what I've read I need to take the cross correlation of these signals and this should give me a peak relating to the time delay. However when I take the cross correlation of these signals I get a peak at 2000 which I know is not correct. What should I do to these signals before I cross correlate them? Just looking for some direction.

EDIT: after removing the DC offset this is the result I am now getting:

Is there a way to clean this up to get a more defined time delay?

EDIT 2: Here are the files:
http://dl.dropbox.com/u/10147354/scope1col.csv
http://dl.dropbox.com/u/10147354/scope2col.csv


texts:
I am a beginner in digatal communications. I have 2 users and I want use a chirp spread spectrum.

Can we take an up-chirp and a down-chirp as spreading signals for the 2 users? and if yes, why? The down-chirp is $x_{down}(t)=x_{up}(-t)$ with $x_{up}(t)$ the up-chirp.


texts: Butterworth filter transfer function in time domain

texts:
Preliminaries
I need to know the response of Butterworth filter in time domain of the following input (at least numerically):
$$
x(t) = \theta(t)(1-\exp\{-at\})\exp\{-bt\}, \quad a,b>0,\quad a\gg b
$$
In order to find the response I decided to find the transfer of Butterworth filter in time domain and then convolve it with $x(t)$ numerically.
But the result I've got doesn't seem to be correct. Here how it went.
Transfer function in $s$-domain
Transfer function of Butterworth filter is the following:
$$
H(s) = \prod_{k=1}^{n}\frac{1}{(s-s_{k})}, \quad s_{k} = \omega_{c}e^{i\frac{2k+n-1}{2n}\pi}
$$
So in order to find transfer function in time domain I used this equality:
$$
h(t) = \sum_{s_{k}}\mathrm{res}H(s_k)\exp\{s_k t\}
$$
Transfer function in $t$-domain
First of all, I found that $\Im{h(t)} = 0$ i.e. it is real. I used this fact as a cross check. I wrote the code that calculates $h(t)$ using complex numbers. I found the following result:

The question is : could somebody verify if it looks right? I haven't found any pictures of the Butterworth transfer function in time domain.
P.S
I could post the code I used to produce those results, but keep in mind it is based on ROOT CERN. You can see the code and tests here.
Update
I have performed one more test: calculation of the magnitude of the frequency response of Butterworth filter of the order $n=3$ and cut frequency $\omega_{c}=2$ rad/s. The result is perfect:


texts: Understanding Voss-McCartney pink noise generation algorithm

texts:
I'm implementing the Voss-McCartney pink noise generation algorithm.

If you follow the link above, you can read:

  from James McCartney  2 Sep 1999 21:00:30 -0600:

    The top end of the spectrum wasn't as good.  The cascade of sin(x)/x
    shapes that I predicted in my other post was quite obvious.
    Ripple was only about 2dB up to Fs/8 and 4dB up to Fs/5.  The response
    was about 5dB down at Fs/4 (one of the sin(x)/x nulls), and there was
    a deep null at Fs/2.
    (These figures are a bit rough.  More averaging would have helped.)

  You can improve the top octave somewhat by adding a white noise
  generator at the same amplitude as the others. Which fills in the
  diagram as follows:

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
x x x x x x x x x x x x x x x x
 x   x   x   x   x   x   x   x
   x       x       x       x
       x               x
               x

  It'll still be bumpy up there, but the nulls won't be as deep.

If I understand it well, this algorithm generates pink noise by adding random (white?) noise sources at different frequencies1

However, I don't fully understand the explanation given in the quote above for the extra white noise generator on the "top row". Can someone clarify how/why it improves the algorithm? Does that make it a good algorithm for pink noise generation for audio applications? Especially, shouldn't I discard the first samples until all the "rows" were mixed into the signal (in the ASCII art quoted above, that would mean discarding 15 first samples)?

1 I'm not sure of the wording here. Do not hesitate to correct me if I'm wrong


texts: Frequency Representation of Downsampled Signal

texts:
I'm trying to follow the steps from Oppenheim for the derivation of the frequency representation of a signal which is to be down sampled by a factor of M.

It makes sense to me that:

But then the book says to make a change of variables from $r$ to $i$ and $k$ where $i$ = $0, 1, ... M-1$, and $k$ = (-$\infty,\infty$). However, I don't really understand how this change of variables comes logically. After making this change of variables it is then easy to see how the downsampled frequency representation is just a sum of scaled and shifted representations of the original sampled sequence. But I would like to be able to make the bridge better between steps and how such a change of variable was concluded.


texts: Multiplication in the wavelet domain, what does it look like in real space?

texts:
Say you wavelet transform a signal, multiply by some values, and untransform. Is there anything we can say about the effect in real space? I know that it's not exactly a convolution, because the convolution theorem doesn't really exist for the wavelet transform. So if it's not a convolution, what is it?

Also, is the effect different if you multiply in an orthogonal wavelet basis (at, say, dyadic frequencies) versus a continuous transform at all frequencies?


texts: Proof that first difference filter amplifies noise

texts:
I'm a bit befuddled by noise's effect on derivative filters.

I've always 'known' that straightforward first difference derivative filters of discrete signals amplifies noise, but I'm struggling to prove that fact.

I thought it was caused by a greater gain for higher frequencies, but when I take the z-transform for a first difference filter I get
$$ H(Z) = 1 - Z^{-1} $$

and if I restrict the magnitude of $Z$ to 1 to get the DTFT and take $|H(Z)|$ the result is
$$ 2 |\sin(\frac{\omega}{2})|$$

which isn't increasing or constant in frequency like a high pass filter.

Could you explain why the noise increases if not for this effect?


texts: Long burst of high amplitude noise removal

texts:
I have a setup with an embedded DSPIC. This is used to link multiple audio inputs to the outputs and changed between them either when told to, or when one is active. This changing takes place on a hardware that is purchased and from what I can tell uses a relay. I have DSP software using a heterodyne filter to remove static noise as the system runs.

However, when the output switches from silent to one output on, there is a large spike in amplitude that lasts for about 3 to 4 milliseconds. This is similar to the click a walkie talkie makes.

I was wondering if there was any way to determine when this happened and then do some sort of amplitude reduction. My sample rate is only 8khz and the signal has been checked on the scope to run at 1khz during this period. I tried doing a 5 point median filter which kind of worked, but on a larger speaker system would still be loud. I also looked into doing a differentiation filter and seeing if the derivative was large. If it was I would cut the next sample point in half and continue. This had no apparent effect.

Any ideas?


texts:
Have following questions about Fixed point conversion format.

Suppose I want to change Fixed Point format as

X1(m1,n1) (→) X2(m2,n2),m1+n1=m2+n2[Word length same],  n2>n1[ Higher Fractionlength]

Suppose  X1=00001100 and suppose I do X1<<(2^3)

Q1)  If I interpret X1 as pure Signed integer, then X1= 12 and above operation means
 Multiply by 2^(n2-n1)=8 so that X2= 01100000 = 96. (Decimal point is irrelevant here?). Is the interpretation correct?
  or
 Assuming Fixed point  notation, above means
 X1(m1,n1) (→) X2(m1-K,n1+k),m1+n1=m2+n2[Word length same],  n2>n1[ Higher Fractionlength]  and k=n2-n1
  For given example, this means X1(3,5)=000.011 00(0.375)(→) X2(0,8)=  .01100000=0.375
  So, the value remained same.

Q2) Can we always interpret this as decimal point being fixed and bits moved left?Is this a logical shift or arithmatic shift?

Q3) Ideally, we would like to preserve resolution and range covered. Under what condition would this be guaranteed for the above shift? for example, consider following conversions:

Consider  X1(3,5)=101.011 00(5.375)(→) X2(0,8) is 101.01100<<3 = .01100000 =0.375.So, we have error.
   Q4) Is this rounding or truncation and why?

Consider  X1(m1n1) = 000...0 . abc......z  , [m1 = all zeroes],  X1 << k = X1(m1-k,m1+k) where k<=m1
            would  not change the value, as bits to the left and right  of decimal point remain same

Q5) So, if fixed point fractional notation is assumed( Positive fraction, with magnitude <1), Left shifts smaller than integer word length would preserve the fractional value?

Suppose I have fixed point number X1(m1,n1) with Word length m1+n1. Suppose  I want to express
it in X2(m2,n2) format so that Word length is reduced ,Fraction length is increased, Integer Word length reduced    but the value remains the same.For example:
X1(3,5) = 000.01100(0.375) -> X2(1,6)
Steps for above could be:

Left shift by n2-n1 bits so that 000 011 00 <<1 = 000 110 00

Keep Least Significant m2+n2 bits so that X2 = 0.011000(0.375) or Keep Most Significant m2+n2 bits so that X2 = 0.001100(0.1875)

Q6) If I was writing a routine that would automatically do conversion, )(say this was Sum result that needed to be converted to output (m2.n2) format automatically), how are what decides whether I keep LS or MS bits or doing a left or right shift?
Suppose I have fixed point number X1(m1,n1) with Word length m1+n1. Suppose  I want to express  it in X2(m2,n2) format so that Word length is increased ,Fraction length is increased but the value remains the same.For example:
 X1(3,5) = 100.01100(-3.625) -> X2(3,13)
and  X1(3.5) = 100.01100(-3.635) -> X2(1,11)

Q7) How can I go about that?  What would be steps if result required increase in Word length and Integer Word Length?
Suppose I need to convert Fixed Point number as per
X1(m+1,n) (→) X2(m,n)[ Word length reduced, Fraction length same, integer Word length reduced]

One of the references I found says do following :
  (x1)&(0xff..e) and keep MS m2+n2 bits
  With X1(1,4) = 0.1011(0.875) this becomes 0.101 = 0.75

Other says: Do left shift by 1 and Keep Ms m2+n2 bits
With X1(1,4) = 1.011(0.875) this becomes 1.011(0.625)

Obviously, i don't perverse the value.
Q8) Is there a correct way to do it and the logic behind it?

I know these are lot of questions but it will help me clear lot of confusing stuff.
Thanks
sedy


texts:
Papers use wireless device-dependent radio-metrics as fingerprints. A radio-metric is a component of radio signal, amplitude, frequency and bandwidth. Each device creates a unique set of radio-metrics in its emitted signal due to hardware variability during the production of the antennas, power amplifiers, ADC and DAC circuits.As a result, radio-metrics cannot be altered post-production, and thus provide a reliable means for distinguishing wireless devices.

In fact, the transmitted RF signal from a wireless device experiences hardware impairments, channel characteristics, and noise at the receiver. Distortions caused by channel-specific and noise-related effects are likely to have a more random structure. Distortions in a feature that are caused by transmitter hardware impairments should manifest themselves consistently across multiple frames from the same transmitter.

device-dependent fingerprints are Modulation based features like carrier Frequency offset, phase shift offset, I/Q origin offset in the constellation plane.

How can I simulate carrier frequency offset in MATLAB?
In the following reference 1, the carrier frequency difference (CFD) is used as a feature for distinguishing wireless devices and it is assumed as Gaussian distribution. Why can we assume that this feature is Gaussian distributed?
Why do noise and the channel not affect this feature?

On Identifying Primary User Emulation Attacks in Cognitive Radio Systems Using Nonparametric Bayesian Classification

//Edit

CFO of different devices in the feature space

//reference

Device fingerprinting to enhance wireless security using nonparametric Bayesian method


texts: Matched filter and cross correlation in digital processing

texts:
I am very new to signal processing, and a bit confused with the matched filter.
Assume that I have a time series and a specific waveform I need to identify in it. By definition:

  a matched filter is obtained by correlating a known signal, or
  template, with an unknown signal to detect the presence of the
  template in the unknown signal.

So I take the window from the sequence, generate the test signal with the same sampling rate, cross-correlate them, find the maximum and threshold it.
Is cross-correlating my measured time sequence with the generated waveform a matched filter (for example, using xcorr in matlab), or there is a specific procedure I must follow to produce a "more correct" matched filter.
Thank you.


texts: shannon interpolation in image processing

texts:
I tried to implement Shannon interpolation on a 2D array.

First, implemented it on a 1D signal, just for sanity-check:

t = -7:7;
t2= -7:0.01:7;
figure
plot(t, sin(pi*t / 3))
upSamp = zeros(size(t2));
for iT = 1:length(t2)
    upSamp(iT) = sum(sinc(t2(iT) - t) .* sin(pi*t / 3));
end
hold
plot(t2, upSamp)

Looks ok.

Then, implementation on a 2D array. The reference array is:

figure
surf(sin(pi*T1 / 3).*sin(pi*T2 / 3))

Then the interpolation is done by a convolution with a Sombrero function:

[T12, T22] = meshgrid(t2, t2);
val = zeros(size(T12));
for iInterp=1:numel(T12)
    xdist = (T12(iInterp) - T1);
    ydist = (T22(iInterp) - T2);
    val(iInterp) = sum(sinc(sqrt(xdist.^2 + ydist.^2)) .* sin(pi*T1 / 3).*sin(pi*T2 / 3), 'all');
end
figure;surf(val, 'EdgeColor', 'none');

seems to work good only on the edges :(

Any suggestion?


texts: Who can tell you how to draw the amplitude-frequency spectrum and the phase-frequency spectrum of a coherent packet of radio pulses?


texts: STFT(Short-Time Fourier Transform) & DFT(or FFT) in Matlab

texts: I want to try STFT & FFT using Matlab. What I wonder is STFT of signal computes the result that FFT(DFT) of each windowed signal and I can see the change of each frequency value over time. If I calculate the average of each frequency over the total time, can I get the same amplitude result with the result of the FFT(DFT) of the whole signal? If not, how can I get the same amplitude result of each frequency for both STFT and FFT(DFT)?


texts:
I am currently implementing acoustic FSK modulation and demodulation. I am not a signal processing guy so any help about the bit timing recovery would be very appreciated. Currently I  implemented the demodulator using two matched filters for each tone (with a difference of phase of $90^{\circ}$ for non coherent detection). Basically the output of each filter are peaks with different amplitudes. I have two questions:

How can I perform bit timing recovery?
What do you recommend for packet synchronization preamble (chirp,barker code, golden code etc)?

I googled out the Gardner algorithm  but I am not sure if it is applicable and how. Consider I am currently working with 2 frequencies with bitrate of 800 bps using the soundcard.


texts:
For a school project, I'm working on a software-defined radio transmitter intended for the HF amateur radio bands. I'm planning to support SSB transmission with the formula

$$f(t) = m(t) \cos(2\pi f_\text{carrier}t) \pm \hat m(t)\sin(2\pi f_\text {carrier}t)$$

where $f_\text{carrier}$ is the IF carrier and $\hat m$ is the Hilbert transform of $m(t)$. Of course, I will use a discrete form of this equation for my implementation, which means I need to calculate the discrete Hilbert transform of $m(t)$.

I've looked around online but can't seem to find a decent C/C++ implementation of the discrete Hilbert transform that would be suitable for my purposes. However, there are plenty of libraries to calculate FFTs.

From my understanding, a discrete Hilbert transform can be calculated by taking the FFT of the signal and multiplying by j to achieve the 90° shift. It suffers from Gibbs' phenomenon, it seems, and might need a wide bandpass filter.

Can anyone tell me if my understanding is correct (or of a good discrete Hilbert transform function)?


texts: Effect of time duration in bandpass filtering

texts:
I am new in DSP, but I am trying to process EEG data to decompose into the bands, theta, alpha, and others.
The raw data is processed:

Remove the DC offset by subtracting the average value from the entire data channels.
High pass filtering.
Both such as in.

The using the bandpass function of Matlab it is bandpassed.
delta = [0.1 4];
theta = [4 8];
alpha = [8 12];
beta = [12 30];
gamma = [30 45];
bands = [delta; theta; alpha; beta; gamma];
fs = 128; % sampling freq of the EEG

% choosing a segment of data

time = 2 ; % secs
fs = 128;
size_t = fs * time;
sample_filt = af3(1:size_t); %  AF3, an EEG node

% band passing
eeg_bands = {};
for i = 1  : size(bands, 1)
    eeg_bands{i} = bandpass( sample_filt, bands(i,:), fs);
end

% code to plot

I performed bandpass for samples of different time duration, however, I saw the time duration changes the amplitudes of bands at the start and end in gamma and beta (high frequency bands).
here, I plot three figures. Figure 1 of long 5 secs (for me, visually it is reasonable). Figure 2, of 2 secs of duration and notice the start of signals beta and gamma are noisy.

Taking in account EEG, I ask you,

is there any ideal time duration (min, or max)  of the signal to make band passing ?
Mi data has a duration of 22 minutes at 128 Hz of sampling, can I performed band pass of all the length or by small segments or windows?

I thank you in advance, please feel free to comment, I need to be right to extract features in the next step.


texts:
I'm trying to implement a simple FHT on .NET platform myself. I follow this document: Optimized fast Hartley transform for the MC68000 with applications in image processing.
In the page 20th, equation (20):

  H(k) = He(k) + [Ho(k)cos(2πk/N) + Ho(N/2-k)sin(2πk/N)] H(k+N/2) = He(k) -
  [Ho(k)cos(2πk/N) + Ho(N/2-k)sin(2πk/N)]

  He(k) is the N/2 point DHT of the even indexed elements of
  H(k) and
  Ho(k) is the N/2 point DHT of the odd indexed
  elements of H(k)

$k$ is zero index, so where is the item $N/2$ if $N/2$ is the length? I know I must have misunderstood something, but I realy don't know what exactly the Ho and He are. I have succeeded implement the FFT equation (in the page 17th). But the FHT, I can not.
Can somebody please help me explain the above equation details?


texts: Determining if a system involving integrals is time varying or time-invariant

texts:
I've been studying for a signals and systems class coming this fall and am trying to understand the following problem.
$$$$
Show whether the following system is time invariant or time varying:
$$S_1:y = T_1[x],\;\; y(t) =  \int_0^t{[x(\tau)]^3\,d\tau},\;\;t\geq 0$$
solution:
$$z(x)=t[x(t-\sigma)] =  \int_0^t{[x(\tau-\sigma)]^3\,d\tau}$$
$$=  \int_{-\sigma}^{t-\sigma}{[x(\tau')]^3\,d\tau'}$$
Where the change of variable $\tau-\sigma \rightarrow \tau'$ is applied. and
$$y(t-\sigma) =  \int_0^{t-\sigma}{[x(\tau)]^3\,d\tau}$$
So I know that the answer is that $z(x) \neq y(t-\sigma)$ and thus the system is time varying. So my question is why doesn't 0 get subtracted by $\sigma$ in integral expressed in $y(t-\sigma)$?

Furthermore if I switch the limits of integration like so:
$$S_2:y = T_2[x],\;\; y(t) =  \int_{t}^{t+1}{[x(\tau)]^3\,d\tau},\;\; t\geq 0$$
Would $S_2$ be time invariant?


texts: How can spectrum of DTFT of sampled band-limited analog signal be in $[-\pi,\pi]$

texts:
Say the analog signal $x(t)$ and its spectrum $X_a(\Omega)$. After sampling with frequency $F_s$ or sampling period $T$ we get

$$x[n] \triangleq x(nT) = x\left(\frac{n}{F_s}\right)$$

and its spectrum $X_d(\omega)$.

The relationship between the two spectrums is

$$X_d(\omega) =\frac1T \sum ^{\infty}_{n=-\infty}X_a\left(\frac{\omega + 2\pi n}{T}\right)$$
$$\omega=T\Omega$$

For band-limited analog signal having spectrum in range $[-B,B]$. $X_a(\Omega)=0$ for $|\Omega| \geq B$ if $F_s$ is above Nyquist rate then
$$X_d(\omega) =\frac1T X_a(\Omega)$$

For digital signal, its frequency is restrict in $[-\pi,\pi]$ but not the case of analog signal. The equation above means that $X_d(\omega)$ is not restricted in $[-\pi,\pi]$


texts: How to zoom an image using DCT transform

texts:
I am familiar with zoom using fourrier transform on an Image.
I am also aware that DCT is used mainly to compress images.
Meanwhile, how can someone use this transform to zoom an image?


texts: How to deal with Imaginary Value of Fast Discrete Curvelet Coefficient?

texts:
I want to classify the mean, energy, std, and entropy of discrete curvelet coefficients(wrapping) with Neural Net. But Neural Net can't deal with imaginary value. I have 2 option, either make it real (remove the imaginary, keep the real) or absolute it ( the imaginary become real). Which one i should choose and why?

e.g. :

abs(3+ 2i) = 5

real(3 + 2i) = 3


texts: Complex Spectral Phase Evolution (CSPE) Performance depending on signal windowing?

texts:
I am look into CSPE. "Signal Analysis Using the Complex Spectral Phase Evolution (CSPE) Method"

The method is simple. It compares the original signal's FFT and shifted signal FFT in phase domain so that it can get an estimate of frequency. The original purpose of this paper is to improve the accuracy. However, I am wondering if it can be used to detect if there is a tone around certain FFT bin.

One way to do this is to get the $\delta$ value for each FFT bin. If $|\delta| < .5 $ indicates a potential tone around the frequency. One simulation is to do it on pure noise. However, I found that the results depend largely on the window you added to your signal. Here is my code:

#!/usr/bin/env python3

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

nfft = 512
nsamples = 513

noise = np.random.randn(nsamples) + 1j * np.random.randn(nsamples)
noise = np.sqrt(.5) * noise

SNR = 10
noise = noise * 10 ** (-SNR/20)

recv = noise # pure noise

s0 = recv[:nsamples-1]
s1 = recv[1:]

S0 = np.fft.fft(s0 * signal.chebwin(len(s0), at = 80), nfft)
S1 = np.fft.fft(s1 * signal.chebwin(len(s1), at = 80), nfft)

#S0 = np.fft.fft(s0, nfft) # square window
#S1 = np.fft.fft(s1, nfft) # square window

SS = np.conj(S0) * S1
aSS = np.angle(SS)
idx = np.where(aSS < 0)
aSS[idx] = aSS[idx] + 2 * np.pi
cSS = aSS * (nfft/2/np.pi)
bSS = cSS - np.arange(nfft)
print(np.sum(np.abs(bSS) < .5)) # estimation of potential # of tones

Some results:

If chebwin is used, I usually got 300 potential tones, which is bad.
If I used squared window, I usually got 20+, which is not bad
If I reduce chebwin's attenuation, the number also reduces.

I can't figure how this can be related to the window function.

Thanks

Yan


texts: Why is it called continuous-time frequency?

texts:
I'm just wondering about the CTFT. My lecturer refers to capital Omega $\Omega$ in the following as being the continuous-time frequency:

Why is it called continuous-time frequency here but in the DTFT the small omega is called discrete-time frequency? CTFT and DTFT both give continuous signals out.


texts: What Is Meant by Clustering of ICA (Independent Component Analysis) Components?

texts: I'm following the EEGLAB tutorial about STUDY creation and manipulation of the datasets, and came across the concept of clustering of ICA components. I'm not too familiar with processing of ICA components other than its physical meaning which is separating a multivariate signal into additive subcomponents. If someone could explain clustering in as simple way to understand intuitively, that would be great!


texts:
I am trying to implement a Chebyshev waveshaper according to the description on the bottom post here:

Add odd/even harmonics to signal?

It looks similar to the algorithm here:
http://www.juce.com/forum/topic/voice-changer-plugin

The results for any odd order setting sound good. But all even order results just end up in massive distortion. I also tried x * abs(x) instead of x², but the result is similar. How should I handle the even order terms? Or do I have to use a completely different method for producing even harmonics?


texts: I am I am working with 'wav' files in audio signal processing. I am able to find out STFT by using the function " spectrogram(file name, window length). Now I want to apply Inverse Short time fourier transform on the above signal. I tried "invspectrogram". but its showing error as Undefined function for input argument type of double. So please tell me how to proceed.


texts: Sampling bandpass signal with spectral replications

texts:
I want to clear my thinking on how to approach problems of the type characterized below.

I guess what we want to do is select a sample rate of a sufficient size such that this signal can be fully captured  therein but also small enough so as not to waste all the space from 0 to, 156? (77*2).

I've seen in some examples that if we select a sampling rate that is some integer multiple of this rate then we can use one that is much lower than the nyquist frequency, taking advantage of the property of replication.

So maybe for part (a), this would just have the same thing repeated over and over every 56 MHz? and How about specral inversion, would every other sample be inverted?

For part (b) maybe it's 56/2 = 28? or maybe it's 77/2?

For part (c) are we meant to think about the

fs' < 2fc-B/m

or maybe the

2fc-B/m > fs > 2fc+B/m+1

I guess m in this situation is 56? And B is 14? But how about the transition of B, it doesn't drop down straight but it flares out at the bottom, how to think about that? And I guess fc = 70?

EDIT: (reaction to Matt L.'s answer)

Is the image below what you had in mind for part (a)?

But why don't we take the bandwidth into account? Furthermore, in the image it seems that the 14MHz indicated there as the bandwidth doesnt completely describe the size of B since it seems to flare out at the ends, how about that?

So can we generalize that equation to be that for any time we are performing that undersampling we can get the location of all the replications by:

'Highest frequency component for which we must account' - ('number of K term' * 'sampling rate')

and sampling rate must be > 2B

Why did you only calculate for two terms? could be go on? can we do k = 3, 4, 5, 6... ?

I guess you mean the answer to part (b) is 14 because it's the smallest positive spectral component right?

& for (c) you mean f1 = 14 MHz = 70MHz - fs


texts: Get the distinctive deltas of quite noisy and very variable / dynamic measurement values

texts:
for my thesis I need to point out the distinctive deltas of my measurements for further use.

The jumps can occur with different heights over different time slots. And the absolute base values can also slowly but not only marginally increase or decrease between ups and downs, caused by other means.

I would love to be able extract these values on the fly on the ARM Cortex M0 that gathers them.. of course delayed, so there are enough measurement values in the pipe to denoise them/calculate medians.. Matlab could also be used (even though I have no experiences in it)..

I thought of something like: Save line when delta is higher than threshold and then sum up everything in between till another delta was undercut (multiple times). Maybe I should also add the "multiple times" argument to the higher than threshold to catch these jumps that consists of multiple small step?

Best, just have a look at this picture:

At the moment I am playing around with the following and get ok values, but it could be better..:

delta > small_threshold? -> start summation. deltasummation > large_threshold after at least x additions? -> wait till delta was < small_threshold for y times -> Result. Its like a low pass + FIR I guess?

I played around with the "small_threshold", "large_threshold",  "x", "y" values till I got the overall best result.. But I think there is room for improvements...

Example Data.csv


texts:
I have got three signals, acc_x, acc_y, acc_z. I need to compose them in order to compute the magnitude of the acc. I would like to ask the following question:

what is better opt1 or opt2? or opt3?:

**opt1**: first apply a filter to reduce the noise to each single signal and then compute the magnitude or,

**opt2**: first compute the magnitude and then apply a filter to the composed to reduce the noise?

**opt3**: Any other good suggestion?

Any of your suggestion would be very appreciated.

Thank you in advance.

Best Regards,
Carlo


texts: Sampling: question regarding sampling period

texts:
what exactly happens when I increase the sampling period ?

Providing some Matlab plots so that I can be more clear.

Using the signal $x(t) = \sin(10\pi t)$ and periods $T_s=0.02 ,\; 0.05 ,\; 0.1$

What happens exactly to the signal when the $T_s = 0.1s$ ?


texts: How to compute the energy of a NON-STATIONARY (transient) random discrete-time signal

texts:
When computing the energy of a NON-STATIONARY (transient) random discrete-time signal $x(n)$, does it make more sense to compute the energy as

$ E=\sum_1^N{x^2(n)}$ over all the $N$ samples

or does it make more sense to compute the energy as the sum of the auto-correlation function values $R_{xx}(k) = \sum_n {x(n+k) \, x(n)}$ , which should be

$ E = \sum_k R_{xx}(k)$ ?

Please let me know what you think about it.

Regards,

E.


texts: Filtering reject two bands (multiband bandstop FIR filter)

texts:
I know the existence of Remez Parks-McClellan, but I want to exploring results with Math....

I have five frequencies on input (check blue plot):

233.08Hz, 246.94Hz, 261.63Hz, 277.18Hz, 293.66Hz

I want to desing a filter that stops two separates bands. I want to remove the second frequency (with bandstop) and remove the fourth frequency (with the same bandstop filter)

But, I want a single pass (apply convolution only one time) to the input!

This is a similar question, but is for bandpass filters:

Using Typical Hamming Window. I have this:

Sampling Rate = $8000\textrm{ Hz}$
Low bands (for 2nd Frequency)

$f_{{Lo}_{c}} = 239.912\textrm{ Hz}$: Frequency $1^{\rm st}$ band cut
$f_{{Lo}_{p}} = 254.178\textrm{ Hz}$: Frequency $1^{\rm st}$ band pass

HighBand (for 4th Frequency)

$f_{{Hi}_{c}} = 269.292\textrm{ Hz}$: Frequency $2^{\rm nd}$ band cut
$f_{{Hi}_{p}} = 285.305\textrm{ Hz}$: Frequency $2^{\rm nd}$ band pass

The calculated size was $N=1961$ (more restrictive size).

I calculate the Coefficients (Java), getBandStopCoefficients is working perfectly (I'm not providing the source code)!!!

http://www.labbookpages.co.uk/audio/firWindowing.html#high

double[] dBSC_B1 = getBandStopCoefficients(N, "Hamming", 239.912, 254.178);  // Band One for Second Frequency
double[] dBSC_B2 = getBandStopCoefficients(N, "Hamming", 269.292, 285.305);  // Band Two for Fourth Frequency

The Response for Second Frequency (Red Plot)

The Response for Fourth Frequency (Green Plot)

I don't want to do something like (but I want the same result):

double[] x;// : input
...
double[] y = convolution(x, dBSC_B1);
y = convolution(y, dBSC_B2);

I was testing something...

when is needed to perform FFT to multiply the data by window

double[] dBSC_Both = new double[N];
for (int i = 0; i < SizeTest; i++) {
  dBSC_Both[i] = dBSC_B1[i]*dBSC_B2[i];
}

//Other manner... Get two BandPass (adding the output, and is subtracted from the input)

double[] dBPC_B1 = getBandPassCoefficients(N, "Hamming", 239.912, 254.178);  // Band One
double[] dBPC_B2 = getBandPassCoefficients(N, "Hamming", 269.292, 285.305);  // Band Two

double[] dBSC_Both = new double[N];
for (int i = 0; i < N; i++) {
  dBSC_Both[i] = 1.0 - (dBPC_B1[i] + dBPC_B2[i]);
}
double[] y = convolution(x, dBSC_Both);
// the BandPass work fairly well, but like double BandStop so bad!

I was trying this ("preconvolution" Check Purple Plot)

double[] dBSC_Both = convolution(dBSC_B2, dBSC_B1);  //
double[] y = convolution(x, dBSC_Both);
// the result was similar poorly than expected (filter effectively attenuated in the bands,
//but not enough, and amplitude permitted in the bands was attenuated)

I was thinking the theory of calculating the coefficients for windows (changing the formula of the sine)

Some clue on how to do this?


texts: ANSI S1.11.2004 Octave-Band Filter implementation

texts:
I am trying to implement some code to firstly generate an Octave-Band and Fractional-Octave Band filter.  I have used MATLAB's octaveFilter method to generate compliant Class-1 and Class-0, 1/1 and 1/3 Octave band filters, but I want to understand how to design these filters from scratch so I can better understand their construction.  Then I want to go on to implement in code.  So the coding environment doesn't matter, but I want to understand the filter construction from the filter definition/algorithm.
My problem is that the documentation for MATLAB's octaveFilter method doesn't explain the algorithm used to create the filters, nor does it explain how the filter is applied to data programmatically.
Can someone please point me to some material that explains how to create these filters, for any centre frequency I choose, and then implement in code?  Many thanks!


texts:
I have this state-space model:

Is this state space model nonlinear? If it is, why is that?


texts: I am at the moment self studying Z-Transform. I need some references on this subject with some applications ( electrical circuits or else). What are good references on this subject. Thanks in advance.


texts: Matlab - Bode plot of Lag Filter + Integrator

texts:
I am rather new to Matlab and I just cant make sense of what I see in the bode plot of the continuous and discrete version of the same function. The bode plot of the continuous function looks as expected. However the bode plot of the discrete version has a phase offset of +90 degrees and the gain stays the same at lower frequencies. It's basically a lag compensator with an integrator. This is not the final result I am going for, but the easiest example I could think of to make the problem as clear as possible. Please note that I am using a rather high sampling rate of ~12 MHz, which is what I need in the final hardware design. The high sampling rate seems to be part of the problem. With lower sampling rates the frequency at which the bode plots differ gets shifted to the left.

Can anyone explain this behaviour and maybe how I can achieve a result that is closer to the continuous bode plot using the high sampling rate?

Code of the example:

FS = 48000*256;
T_sample = 1/FS;

accu = tf([0 1], [1 0]);
lag = tf([10 1], [100 1]);

lag_d = c2d(lag, T_sample, "zoh");
accu_d = c2d(accu, T_sample, "zoh");

bode(lag*accu, lag_d*accu_d);


texts:
I am trying to use a FIR filter making software to create a moving average filter and I'm a little confused on what I should put in the coefficients part? I read that the coefficients are always 1 for moving average so this is what I did - is that correct? (Picture below)
Other than the coefficients, I just wanted to double check (because I am very new to working with filters):
Moving average filters are symmetrical, have no L-th band filter, and are single rate filters right?


texts: Equation for half-damped triangle signal

texts:
I have triangle signal, which is damped only at the beginning of ramping phase (blue line in figure).

I want to fit this type of signal. Just now I use symmetric equation, which gives me a damp at the end of ramping also:

$$ y=C + \frac{A}{\pi}*\sum_{n=1}^{n=k}{\left[ {-1^n}\frac{\sin(2\pi nf(x-x_0))}{n}\right]} $$

Could You help me to find an equation which will give damping only from one side? Do you know similar signals? I will be good for me to read how to deal with them.

To clarify:

I have a set of signals, which looks as blue line.
I need to fit them. Signals has "damping" peak only at the beginning of ramping('wanted peak'). Solution I have, gives me also an "unwanted peak".
I am searching for better solution. Which will gives me a green line, now i draw green by hand.


texts:
If we want to change resolution of an  image,we can achieve that by increasing or decreasing the number of pixels

Is there any other way ? Can we change dimensions of an indvidual pixel,rather then changing the number/count of total available pixels?

by changing size of pixel,i mean changing the width and height of array/matrix containg pixel value


texts:
I was unable to solve a convolution question.Question is attached herewith.

I don't know how to initiate for solving the problem and what is the final expression that prove the convolution of  a top-hat function with itself is the triangle function.


texts:
in wireless communication,compared to the OFDM system , why does the SC-FDMA (single carrier FDMA) system has a lower PAPR (peak to average power ratio)?

regards,
phani tej


texts:
I have this binary image with some blobs of random sizes.

I would like to convert the blobs into circles of uniform diameter. it doesn't really matter what diameter but they need to be circles and uniform.
Here is what I am looking for.

I am working with tiff sequences in ImageJ. I tried erosion and dilation. I tried Ultimate point to reduce to a point and then dilate up but that creates a square. I tried distance map but didn't get very far. Although i am pretty sure distance map is the way to go. None of my current approaches are satisfactory.

If I could do

   #pseudo-code
   if area(blob) > x
   then make circle(diameter=y)
   else remove blob

would be even better.

If anyone has a solution using ImageJ or Python SciPy or even R, please let me know. I would very much appreciate that. Thanks.


texts:
General Question

Given an audio file, say a 16-bit wav, what are some standard methods to preprocess a spectrogram of this wav so that it may be fed into a neural net?

Context

In Lee et al's 2009 NIPS paper on applying deep belief networks (DBNs) to speech recognition, they briefly mention in Sections 2.2 and 3.1 preprocessing a spectrogram using PCA whitening.

NB: I don't have image-processing background as applied to machine learning...

Question

How may we reproduce the training data that is fed into the DBN in the above NIPS article, section 4.1 for speaker identification in the TIMIT dataset?

"We trained the CDBN on unlabeled TIMIT data (as described in Section 3.1)...We computed the first and second-layer CDBN features using the
spectrogram as input."

Section 3.1 mentinons 20ms window sizes with 10ms overlaps to generate the spectrograms and then applying PCA whitenting. Does this spectrogram -> whitening process produce an $n \times  m$ or $n \times m \times t$ array for each speaker? Or do we get an "image" that is somehow fed directly into a (C)DBN?

Outlook

Based on the NIPS paper, such preprocessing seems standard but I'm (obviously) new to using spectrograms in neural networks.

Any python references on spectrograms are appreciated as related to these questions, for example with librosa, pylab, or this tutorial.

After generating a spectrogram, perhaps the general steps in this PCA whitening overview can be applied to produce relevant input to a DBN?

Example

For many supervised classification tasks, we may fit a model via fit(X, y) where $X$ is an $n \times k$ matrix of $n$ samples with $k$ features and $y$ is an $n \times 1$ vector where each row is a corresponding label.

Using librosa's melspectogram method on an audio file, we produce an $m \times t$ matrix. With, say, $n$ labeled audio files we would have $n$ training samples each of shape $m \times t$ resulting in an $n \times m \times t$ matrix $X$ (and $n \times 1$ matrix $y$). How may we feed this $X$ into a DBN or neural network in general?


texts: Digital Signal Processing Filters Problem

texts:
Hello I need some help with this question.
$$[] = _1[ − _1] + _2[− _2] + _1 [− _3]$$
Set $a_1, a_2$, and $b_1$ to any positive or negative real numbers of your choice. Set $k_1, k_2$, and $k_3$ to any positive integer numbers of your choice ($k_1$ can be zero). Present this equation with the numbers you have chosen. Present the transfer function $()$ of your filter
I have chosen the values so my equation looks like this
$$y[] = 0.4[] + 0.6[− 1] + 0.8[− 2]$$
I have worked out $H(z)$ to be $$\frac{0.4 + 0.6z^{-1}}{1-0.8z^{-2}}$$   is this correct?
My questions are:

Is this filter an FIR or IIR type?
How do you determine the poles and zero?
How do you determine if the filter is stable or unstable?

Thanks!


texts: Maintain constant total amplitude of sum of sines

texts:
Given a set of N sine waves, each with their own amplitude and frequency, is there a straightforward way to select their phases such that the sum over all sines maintains a constant amplitude, given that the frequencies are spaced equally apart by df?

If all phases are 0, then at t=0 all the sines will constructively interfere to make a huge spike, then they will destructively interfere until they become commensurate again at t~1/df. My goal is to select a set of phases for the different sine waves which will leave the resultant sum with roughly constant peak-to-peak amplitude from t=0 to t=1/df. Picking all random phases does do a pretty good job, but ideally I'd like to find a way to find an optimal set of phases, or at the very least come up with a metric to quantify "how good" a given set of phase is at accomplishing my goal, which I can then optimize on.

edit: To clarify I'm not looking to have the sum of sines equal a constant value, but rather have the resulting periodic function have a roughly constant peak-to-peak amplitude. For instance, for a given set of frequencies all with equal amplitudes, if all phases are set to 0 then the resulting sum will look like:

If the phases are set uniformly randomly in the range [0 2pi] then the result is:
.

Clearly, in the random-phase case the peak-to-peak amplitude of the resulting oscillatory function is roughly constant over the given time series, while in the 0 phase case it spikes often at times where commensurate sines add constructively. I'm looking for a way to find an optimal set of phases which keeps the peak-to-peak amplitude roughly constant, without just generating a million different random combinations and picking the best one.

Any help would be greatly appreciated.


texts: How to reconcile "bandwidth" in the Shannon-Hartley Theorem with a spread-spectrum principle

texts:
I'm trying to understand something about channel "bandwidth" $B$ in the Shannon-Hartley Theorem: $$C = B \log_2 \left( 1 + \dfrac{S}{N} \right) $$

Suppose I have a bitstream encoded as a signal that occupies some bandwidth $B_1$ which is less than channel bandwidth $B$, and I have additive white Gaussian noise $N$. Shannon-Hartley tells me that I'm not using the communication channel efficiently; I should be bumping up the bandwidth of my encoder.

One way of doing that is to use spread-spectrum... but I've read several places (including this paper) that SS doesn't help against wideband noise; it only allows you to increase SNR if you have narrowband noise.

Here are my questions:

If I'm only using bandwidth $B_1$, doesn't that mean I can filter some of the noise out? And if so, would I be able to get $S/N$ down up? Wouldn't that be an alternative to occupying the full channel bandwidth?
Is spread-spectrum a bad idea to try to occupy the full channel bandwidth?


texts: Intuitive definition of ergodicity for random signal

texts: Is it possible to define the ergodicity of the random signal in an intuitive sense without using any statistical reference?


texts: Why are sidebands generated in AM and FM?

texts:
When the signal is modulated onto the carrier in the electromagnetic spectrum,
that signal occupies the small portion of the spectrum surrounding the carrier frequency. It also cause sidebands to be generated at frequencies above and below the carrier frequency.

But how and why are those sidebands generated in AM and FM and why are there so many sidebands generated in FM while just two are generated in AM ? Please provide a practical example, as I already know how they are generated mathematically.

What I know is, in the time domain, when the original signal is put into the carrier signal, it is actually multiplied with the carrier signal which means that in frequency-domain the original signal is convolved with the carrier signal. Those two Sidebands in AM are actually the Fourier transform of the carrier signal.

Is this correct?


texts:
Scheduler schedules the resource for each UE. How often do a LTE eNodeB scheduler schedule? Is the time between each scheduling fixed?
 any idea/reference will be great help on this regard. Thanks.


texts:
I am new in the field of systems and signals, and I have a rather basic for the majority of the group, question:

Can we find the impulse response function of homogeneous ODE, instead of its zero-input response?

for example, we have the following 2nd order homogeneous ODE:

$a_{2}\ddot{x}+a_{1}\dot{x}+a_{0}{x} = 0$

, where the output is the $x(t)$ given the initial conditions

I understand that if it were:
$a_{2}\ddot{x}+a_{1}\dot{x}+a_{0}{x} = f(t)$

the $f(t)$ would be its input, $x(t)$ its output given the input, and we could find the impulse response by replacing the $f(t)$ with $\delta (t)$.

Now that the input is zero, how can we find what the output would look like with respect to any input?

Is this:

$a_{2}\ddot{x}+a_{1}\dot{x}+a_{0}{x} = \delta (t)$

even allowed, for an initially homogeneous equation?


texts: What is the Fourier Transform of an integral

texts:
If we have a function that is an integral over the interval 0 to T

Is is correct to say that its Fourier Transform is :


texts:
I am applying a flat-top window to signal because I am very interested in its actual amplitude in the time-domain. I have read that overlapping will improve results. I am using MATLAB and the following is the code that I am using to compute FFT with flat-top window:

time = 0:0.1:500;
data = 1.8*sin(2*pi*0.2003*time);Fs = 10;
N = length(data);
w = flattopwin(N);
wdata = data(:).*w;
ws = sum(w);
Y = fft(wdata)/ws;
freq = 0:Fs/N:Fs/2-Fs/N;
freq = freq';
amplitude = 2*abs(Y(1:floor(N/2)));
plot(freq,amplitude);

Can someone please help me to modify the code so that I can include overlapping? I am relatively new in this field and am still confused. Thanks.


texts:
Possible Duplicate:
  How is Weighted Thresholded Histogram Equalization different from Gamma Correction Image enhancement

In case of Weighted Thresholded Histogram Equalization we use:-

calculate, the Weighted Probablity Density Function for Image using,

Pwt(k) = ((P(k)-Pl)/(Pu-Pl))^r* Pu,

Where Pu is the highest probablity of a intensity in a Image.
      Pl is the lowest probablity of a intensity in a Image.
      P(k) is the Probablity Density Function for the Image,ie the no of pixels having the   intensity k in a image.

then perform Histogram Equalization on the image using,
calculating Cumulative Distribution Function

Cwt(k) = m=0∑k Pwt(m);
Im(i,j) = Wout*Cwt(F(i,j));

For r<1 the Power Law Transformation function will give a higher weight to the low probablities in the PDF than the high probablities.

How is this procedure conceptually different from the Gamma Correction technique , explain if similiar how much similiarity is there or if something is different , Explain where does it deviate from Gamma Correction method for Single Image Enhancement.

Does this Procedure Passes for Being Called as Gamma Correction in some Aspect OR Do some Changes are Required.


texts: Inverse z transform - contour integration

texts:
Here is my task:

Find inverse z transform of $X(z)=\frac{1}{2-3z}$, if $|z|>\frac{2}{3}$

I need to find it using definition formula, $x(n)=\frac{1}{2\pi j}\oint_{C}^{ } X(z)z^{n-1}dz$. How can I find it?
Result is $-\frac{1}{3}(\frac{2}{3})^{n-1}u(n-1)$


texts: Understanding why spectral flatness cannot be computed

texts:
I have a real signal recorded 512Hz for 1s. After resampling the signal at 256Hz, I would like to compute its spectral flatness. To do so, I start by taking its discrete Fourier transform with fft = scipy.fft.rfft(x) and then I compute the power spectrum with pw = np.abs(fft). When I compute the geometric mean of pw an error is raised because that the last coefficient of fft is equal to 0. I am trying to understand why.

This last coefficient returned by fftpack.rfft is defined as:
$$
\begin{align*} \mathrm{Re}\big( y(n/2)\big) & = {} \mathrm{Re}\Big( \sum_{k=0}^{n-1} x_k e^{-jk\pi}\Big) \\[2mm]
 & = x_0 - x_1 + x_2 - x_3 + \ldots
\end{align*}
$$
where $(x_k)_{0 \leq k \leq n-1}$ is my resampled signal. A quick check that this quantity is close to 0 : np.sum([x[k]*((-1.) ** (k % 2)) for k in range(x.size)]) gives me 5.68e-14, which is not 0 but quite close.

What does that say about my signal ? Usually the spectral flatness is defined without assuming that some coefficients in pw could be zero. So, am I doing something wrong (or not doing something I should do) ?

Edit : is there a convention which consists in saying that the flatness of the signal is 0 whenever a single "coefficient" returned by rfft is zero ?


texts: Vector of random numbers with fixed sampling rate

texts:
How to create a vector of random numbers(in Matlab), when each value is [0,1] .

And the sampling rate should be fixed ?

thanks


texts: What are variance and bias in spectral estimation (specifically periodogram spectral estimation)?

texts:
So far, I have read that all the non-parametric estimation techniques decrease the frequency resolution in order to decrease the variance in the spectral estimate

What is the general "overview" meaning of variance and what is the advantage in decreasing the variance of spectral estimate?

The PSD can be calculated from the autocorrelation function of a signal by taking its Fourier transform, called the periodogram. In general the variance of the estimate $P_{xx}(f)$ does not decay to zero as $N$ tends to infinity.  Thus, "the periodogram is not a consistent estimate of the true power density spectrum", i.e., it does not converge to the true power density spectrum. Thus, the estimated spectrum suffers from the smoothing  effects and the leakage embodied in the Barlett window. The smoothing and leakage ultimately limit our ability to resolve closely spaced spectra. The methods Barlett along with Blackman, Tukey, and Welch are classical methods and make no assumption about how the data were generated and are hence called non-parametric methods.

In this method for reducing the variance in the periodogram involves three steps:

First the $N$-point sequence is subdivided into $K$ non-overlapping segments, where each segment has length $M$.
For each segment we compute the periodogram.
Finally, we average the periodograms for the $K$ segments to obtain the Barlett power spectrum estimate.


texts: Causality and ROC of a stable LTI system

texts:
So I am looking at a stable LTI system whose input is $x[n]$ and output is $y[n]$. The equation relating the two is here:
$$
y[n-1]-\frac{10}{3}y[n]+y[n+1]=x[n]
$$

I was able to compute its system function $H(z)$ to be:
$$
H(z)=\frac{z^{-1}}{z^{-2}-\frac{10}{3}z^{-1}+1}
$$

I then performed partial fraction expansion to end up with:
$$
H(z)=\frac{\frac 38}{\left(1-3z^{-1}\right)} - \frac{\frac 38}{\left(1-\frac 13z^{-1}\right)}
$$

What I am struggling with is I can't decide what its ROC should be since the system is specified to be stable which means it has to include the $j\omega$-axis. I also need to determine if the system is causal or not. Since a system can only be causal and stable when all its poles are on the left hand plane, I am not sure.

Any help determining both the ROC and causality of the system would be very helpful.


texts: Why are we always interested in mean-squared distortion?

texts:
When dealing with quantizers, and in many other communications problems, the interest is usually on the mean-squared distortion or mean-squared error, rather than mean absolute error or anything else. Gallager explains the reason here, probably perfectly if I were to extrapolate from the quality of the rest of his lectures, but I don't get it.

He says that the mean-squared distortion of the quantizer maps nicely to the mean-squared distortion of the waveform we sample&quantize, and this mapping wouldn't be as nice if we used another metric such as mean absolute distortion. Why is that? Shouldn't the mean-absolute difference between the original and the sampled-quantized waveform also map to the mean-absolute difference of the quantizer? What makes mean-square different in this sense?


texts: How can a medical doctor use the information given by someone working in medical image computing(shape analysis)? Details follow

texts:
Posted also here, where it's been put on hold. I did modify it afterwards though.

https://biology.stackexchange.com/questions/43997/how-can-a-medical-doctor-use-the-information-given-by-someone-working-in-medical

I'm an academic researcher in the US working in an area called "statistical shape analysis", which is intimately connected with medical image computing. It deals with how to mathematically compare two images, or predict (regression problem) how the future images (say MRI image of Alzheimer's patients, or images of protein structures) will be.

Let's say I'm working on statistical modeling of hippocampus (a part of the brain) to infer about which stage the Alzheimer's disease is. So I obtain the segmented MRI images, perform statistics on it to predict future segmented images. But I can't tell how exactly the image itself will look like, but only how the segmentation will be like. I hand this prediction to a doctor. How will that be helpful to a doctor who's treating Alzheimer's? Of course, intuitively it's clear it's helpful, but I'm asking what exactly the doctor will do with these kind of computational data.

Also, could you mention some hospitals in Europe where they use informaticians/statisticians to collaborate with doctors for the above kind of work?


texts: Questions in a Character Degradation Model Paper

texts:
I have questions in the ‘Noise region definition’ and ‘Noise generation process’ of the paper “A character degradation model for grayscale ancient document images”.

In Noise region definition, g controls the flatness of the regions. What does it exactly mean? How can we say that a noise region is flatter compared to another noise region?

Below is the illustration of the ellipse noise region within the document image. The green ellipse shape is defined to be the noise region.

According to paragraph below, the average value bj of its 8-neighbours (in the initial grayscale image) is used for calculating values of all pixels in the line CiBj.

Does the average value of bj  calculated by averaging the greyscale value of the adjacent pixels in north, northeast, east, south east, south, south west, west and north west?

Does the formula for getting Pk is below.

Please refer to the paragraph below:

I just want to assess , if my comprehension in reading the article is right. Thanks.


texts:
I am working on OFDM over harsh channel, so the channel coding is an essential thing for achieving a reliable communication.
According to my reading, LDPC is almost the best channel coding we can use for channel coding, but I have a question regarding that kind of coding.
Assume we are using LDCP with rate $1/2$ which means that communication systems uses the half of data rate for coding and the other half as data.
If we have the OFDM symbol have $1024$ sub-carriers, with modulation order $M=4$ = which gives us a total bits of symbol $2048$. In that case, will the LPDC will take the whole $2048$ together? I mean it will take $128$ bits from the whole bits for coding and the rest for data OR it will take group by group, which means for example group of $6$ bits, and then code them $3$ for data and $3$ for parity or coding.
Which one is right ?


texts: Why is this FFT code not working to extract frequencies from a wav audio?

texts:
I'm quite new to audio processing and i'm trying to extract frequencies from a wav file. As research papers stated i put the wav audio into a byte array and segmented and trying to apply FFT to get the frequencies from the audio. However the problem is even when i enter a audio track with only one note played (With no noise) i get a lot of frequencies as a result where as i'm supposed to get one.

My guess is something is wrong with the FFT code because the audio gets segmented properly.

Thanks everyone in advance,

My fft code

public void getFFT(double[] x, int a) {

double amp[] = new double[512];

double firstSum;
double secSum;
double twoPInjk;
double[][] s = new double[512][2];
for (int j = 0; j < 512; j++) {
    firstSum = 0;
    secSum = 0;
    FreqArr[j] = (j * sampleRate) / (512);
    for (int k = 0; k < 512; k++) {
        twoPInjk = ((2 * Math.PI) / 512) * (j * k);
        firstSum += x[k] * Math.cos(twoPInjk);
        secSum += x[k] * Math.sin(twoPInjk);
    }
    amp[j] = Math.abs(Math.sqrt(Math.pow(firstSum, 2) + Math.pow(secSum, 2)));
}
highest(amp, a);
}

public void highest(double[] amp, int a) {
String note;
double all = 0;
double max_magnitude = -1;
int max_index = -1;

for (int i = 0; i < (512); i++) {

    all = all + amp[i];

    if (amp[i] > max_magnitude) {
        max_magnitude = amp[i];
        max_index = i;
    }
}
double freq = 0;
int cnt = 0;

freq = max_index * (sampleRate / 2) / (512);
System.out.println("Maximum_index: " + max_index);
System.out.println("Sample rate: " + sampleRate);

  System.out.println(freq + " is the frequnecy");
}

my result (for the 1 second audio with only one note played)

Sample rate 44100 My data size 70376 My myBitsPerSample 16 My channels 2 My myByteRate 176400 My Length 0.7983219954648526 0.7983219954648526 44100.0 Max_index: 56 Sample rate: 44100.0 2411.71875 is the frequnecy Max_index: 256 Sample rate: 44100.0 11025.0 is the frequnecy Max_index: 256 Sample rate: 44100.0 11025.0 is the frequnecy Max_index: 129 Sample rate: 44100.0 5555.56640625 is the frequnecy Max_index: 511 Sample rate: 44100.0 22006.93359375 is the frequnecy Max_index: 256 Sample rate: 44100.0 11025.0 is the frequnecy Max_index: 152 Sample rate: 44100.0 6546.09375 is the frequnecy Max_index: 422 Sample rate: 44100.0 18174.0234375 is the frequnecy Max_index: 26 Sample rate: 44100.0 1119.7265625 is the frequnecy Max_index: 462 Sample rate: 44100.0

and it goes on.


texts: Please help me understand this FM carrier frequency inequality

texts:
The link of the statement is evidenced in this book:

I don't understand how equation 85 come? I try to use Carson rule to see this inequality but turns out I have $2 \Delta f$ instead of $3 \Delta f$.

So someone could help me to understand this?


texts: Does anyone have a simple algorithm for computing a reasonably accurate arcsine? By "simple" I mean some sort of polynomial that requires <= 5 multiplies per output sample. And by "reasonably accurate" I mean an algo whose error is no more than 10% when the input argument is close to plus or minus one. I searched the web for a while but found nothing immediately useful.