texts: Finding good feedback coefficients for basic reverb

texts:
I'm currently working on some reverb algorithm code for the Gamebody Advance to get a bit of that ambient feeling for the music instead of having plain dull sound. A few challanging restrictions are the low CPU speed of the GBA (~16 MHz ARM CPU) and the 4KiB of memory I have available for the reverb buffer @ 31536 Hz. Perhaps I could do the reverb at half the rate to save some memory. Goal is to keep the CPU budget for this under 5% and I'm well aware of this being a very eager goal.
I've tried a few simple things like a simple feedback line but the results didn't turn out that great and sounded more like echo than actual reverb. However, I could get the ASM of this to work at about 2-3% CPU time.

So before trying to implement anything on hardware I wrote a little C++ application to try out a few algorithms and without the quirks of having no floats on the GBA. This is the core piece of code I've came up with:

// processing is done inplace on "buffer"
// samples ares stored L/R alternating
//
// reverbBuffer is a simple mono buffer
// buffer wraparound is already taken care of
//
// My reverb buffer positions are selected by being 'approximately' coprime
//
// NBPOS = 3 (number of taps on the buffer)

for (size_t i = 0; i < todo; i++) {
    float dryL = buffer[0];
    float dryR = buffer[1];
    float wetL = 0.f;
    float wetR = 0.f;
    for (int i = 0; i < NBPOS; i++) {
        wetL += coeffL[i] * reverbBuffer[bufferPos[i]];
        wetR += coeffR[i] * reverbBuffer[bufferPos[i]];
    }

    *buffer++ = dryL + wetMix * wetL;
    *buffer++ = dryR + wetMix * wetR;

    reverbBuffer[bufferPos[0]] = ((dryL + dryR) * 0.5f) * dryMix + ((wetL + wetR) * 0.5f);
    for (int i = 0; i < NBPOS; i++)
        bufferPos[i]++;
}

For testing I am using it on a few wave files @ 44.1kHz with an appropriately longer reverb buffer.

So, the results we're ... mixed. I've added some logic to create random values for the feedback coefficients (-1.0 to 1.0), the "dryMix" and "wetMix" variables (0.0 to 1.0). Depending on how lucky I either got terrible metallic oscillation or reasonably good results (still nothing I'd call perfect).

Is there any method of finding good coefficients for my feedback taps? Apart from bad sound I also had problems with instability and self oscillation. I don't have any professional DSP knowledge but I'd guess that the sum of the absolute values of all my coefficients should not exceed 1.0 (however this didn't seem to be always true for negative coefficients).

Also, if there is no way to calculate these coefficients with a lot of sense, is there at least a way of randomly generating these in a way they don't cause self oscillation?

PS: I've googled quite a bit about reverb and various things but I didn't find an algorithm which perfectly suits my case. Usually they are too basic (single tap feedback line) or they aren't fast enough (because I need something VERY FAST with low memory footprint).


texts: Why downsample exhibit wider bandwidth when doing discrete time fourier transform?

texts:
Based on Nyquist sampling theorem, sampling would have a convolution in frequency domain, and reasonably the bandwidth of each convolution would be as that of the original signal.
Nevertheless, the discrete time Fourier transform of a signal ,which has been  downsampled , has the $M$ multiples bandwidth of each convolution.($M$ is the downsample factor). I could not have got the idea, what is the true relation between time/frequency domain with underlying discrete time Fourier transform?

Sampling gets convolution in frequency domain:

Downsampling gets more wider bandwidth:


texts:
I have a 6 state Kalman Filter (Unscented).

When I use a diagonal matrix only for Q (i.e only the diagonal has covariances), I get a "smooth" plot of estimate against actual.

If I use the entire Q matrix (i.e include all covariances) I get an estimate that oscillates around actual wildly.

Is it incorrect to just use diagonal?

ADDITIONAL:

I am using 6 correlated Equity Stocks Price timeseries. In constructing the covariance matrix I use the difference between consecutive prices. I do NOT log the data.


texts:
I have samples of a signal that contains an AM signal added to an FM signal with the same carrier frequency $f_c$. The AM signal is a sinusoid with frequency $f_{AM}$ and the FM signal is a sinusoid with frequency $f_{FM}$.

I know how to recover $f_{AM}$ and $f_{FM}$ when I have only an AM signal or only an FM signal, but my techniques are not working when I try to recover $f_{AM}$ and $f_{FM}$ from the sum. For the AM signal, I'm taking the FFT and finding peaks at $f_c$ and $f_c + f_{AM}$. For the FM signal, I'm finding the time-derivative of the signal and taking the FFT, which gives me a peak at $f_{FM}$. When I take the FFT of the sum, I can see peaks at $f_c$ and $f_c + f_{AM}$, but when I take the FFT of the time-derivative of the sum, I don't get a single neat peak anywhere near $f_{FM}$.

What should I do to recover both $f_{AM}$ and $f_{FM}$?

EDIT:
A partial solution I've found is to proceed in two stages. First I find the AM component by taking the FFT (this gives $f_c$ and $f_c + f_{AM}$). Then I eliminate the AM component by dividing each sample by its magnitude. Taking the fft of the time-derivative of this normalized signal gives me a spike at $f_{FM}$, but it also gives me a lot of other spikes. Any suggestions on how to make this give a cleaner answer?


texts: How to find Sampling frequency of a video signal

texts: i took the Digital Image processing subject in college and i am a little confused that how i can find out sampling frequency of a video signal given its bandwidth. For example what will be the sampling frequency of a video signal of bandwidth 18kHz or 15kHz.


texts: Is there any computational method to prove whether a series is stationary or not?

texts: I have a discrete series $x[n]$. It is extracted from real life and I do not have probability distribution of each value $x[n]$. Is there any computational method to prove whether the series is stationary or not?


texts: FFT showing evenly spaced signals when tuned to ~56 MHz down to 25.680

texts:
I'm running GQRX on Mac OS X 10.11.6 hooked up to an RTL-SDR dongle.  I noticed when I'm tuned anywhere between ~56 MHz on down to 25.680 MHz I see evenly spaced signals on the FFT.  But when I tune to 25.679, they disappear completely.

On the higher end, the signals are faint (here I am tuned to 56 MHz):

They become more pronounced as I work my way down (here I am tuned to 43 MHz):

They are most pronounced when I'm tuned to 25.680 MHz.  But when I tune to 25.679 MHz they disappear entirely.  Here I am switching back and forth between 25.679 and 25.680 MHz:

What's causing this?  Since I see signals below 25.3 MHz when I'm tuned to 25.680 but they disappear when 25.679 I tend to think this is a bug.  If so, is it hardware or software related?

EDIT:

The signals are more pronounced tonight.  Each one has a valley surrounded by two peaks:

The one I'm tuned to in the image above has peaks at 27.03 and 27.05 MHz.  The first peak is at about -60 decibels and the second one is at -50.


texts: Get FIR Hilbert coefficients from Matlabs envelope() function

texts:
How can I get the used Hilbert filter coefficients from the Matlab function envelope()?

filter_length = 10;
[up1,lo1] = envelope(input_data, filter_length,'analytic');


texts: Deriving the integration property of the Fourier Transform

texts:
I want to derive the property of the Fourier Transform that states that if $X(j\omega) = \mathcal{F} (x(t))$ then $$\mathcal{F} \left( \int_{-\infty}^{t} x(\tau) \mathrm{d} \tau \right) = \frac{1}{j\omega} X(j\omega) + \pi X(0) \delta(\omega).$$
I started as follows, since $$x(t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} X(j\omega) e^{j\omega t} \mathrm{d}\omega,$$ integrating both sides gives us $$\int_{-\infty}^{t} x(\tau) \mathrm{d} \tau = \int_{-\infty}^{t} \frac{1}{2\pi} \int_{-\infty}^{\infty} X(j\omega) e^{j\omega\tau} \mathrm{d}\omega\mathrm{d} \tau$$ $$\int_{-\infty}^{t} x(\tau) \mathrm{d} \tau = \frac{1}{2\pi} \int_{-\infty}^{\infty} X(j\omega) \left(\int_{-\infty}^{t} e^{i\omega \tau} \mathrm{d}\tau \right) \mathrm{d}\omega.$$

How do I evaluate the inner integral inside the parantheses? My gut feeling says that it's related to the Fourier Transform of the constant function $y(t) = 1$ but I can't relate it to that since the upper limit is $t$ rather than $\infty$. How do I reach the result above?

Thank you in advance.


texts: what is silence in audio .wav file and how to compute decibels?

texts:
I'm a beginner of audio signal processing, and trying to understand some basic things about audio signals.
I primary work with .wav files, so please allow me to use wav as my data to describe my questions.
Q1: Sample values in .wav file represent the waveform amplitude?
When we load data from .wav file, we will have the audio sample values (typically signed ints), and I searched on google, these values are the amplitude of the audio signal, right?
See example code
In [69]: from scipy.io import wavfile

In [70]: frame_rate, values = wavfile.read("data/test.wav")

In [71]: frame_rate
Out[71]: 48000

In [72]: values
Out[72]:
array([[  0,   0],
       [  0,   0],
       [  0,   0],
       ...,
       [385, 247],
       [412, 266],
       [359, 232]], dtype=int16)

Q2: What is silence?
If the answer to Q1 is YES, then given bit-16, I tried to create a numpy array with all values of 10000, expecting it should be loud sound. However, to my surprise it's just no sound at all. Why? I though value=0 is silence, but it's obviously not.
See example code
import numpy as np
from IPython.display import Audio

# 100k frames
x = np.arange(0, 100000)

# sine wave, max value=10000
sine = np.sin(x) * 10000
Audio(sine, rate=50000)  # i.e. 2sec audio, and I can hear the sound

# dummy, all values = 10000
dummy = np.ones(x.shape[0]) * 10000
Audio(dummy, rate=50000)  # same 2sec, but silence

Q3: How to compute decibels from the .wav amplitude values?
Given Q2, now I'm lost and don't know how decibels look like given amplitude values, I thought bigger amplitude values should mean higher decibels, but now it's apparently not the case.


texts: struggling to understand why Fourier basis is orthogonal

texts:
Studying DSP on my own time on Coursera.

Was given a proof to why the Fourier basis is orthogonal, but I can't figure it out. Here is how it is proof goes.

Consider the Fourier basis
$$
\left\lbrace \mathbf w^{(k)} \right\rbrace_{k=0,...,N−1}
$$

defined as:

$$
\mathbf w^{(k)}_n = e^{−j\frac{2\pi}{N} nk}
$$

Let us compute the inner product, that is

\begin{align}
<\mathbf w^{(k)}, \mathbf w^{(h)} > &= \sum_{n=0}^{N-1} \mathbf w^{*(k)}[n] \mathbf w^{(h)}[n]\\
&= \sum_{n=0}^{N-1} e^{j\frac{2\pi}{N} nk} e^{−j\frac{2\pi}{N} nh}\\
& = \sum_{n=0}^{N-1} e^{-j\frac{2\pi}{N} n (h-k)}\\
&= \begin{cases} N & \text{if}\quad k=h\\[2ex] 0 & \text{otherwise} \end{cases}
\end{align}

I do understand all of the derivation except the last step. How does one arrive at the conclusion that when $k\neq h$ the inner product is $0$. Would appreciate any help with that.


texts: Calibrating cameras with different resolutions

texts:
I have two cameras with very different resolutions. One camera is the RGB camera of the microsoft kinect sensor with a resolution of 1920x1080 and the other one is a very specific camera with a resolution of 480x320 (not the IR sensor of the kinect).
Calibrating the intrinsics of the cameras seperately works out just fine. I am getting reprojection errors of around 0.1-0.15 pixels. But when I try a stereo calibration on the sensors, I get very high errors of 1.0 pixels and higher.
For calibration I used OpenCV. First I was trying to use Matlabs stereo calibrator, but it didn't work for cameras with different resolutions.
Do you have any hints on how to improve my results? Is there anything special that I should do different when working with very different cameras?


texts: Difference equation when transfer function expressed as poles and zeros

texts:
The transfer function $H(z)$:

$$
H(z) = \frac{Y(z)}{X(z)} = \frac {b_0 + b_1 z^{-1} + b_2 z^{-2}} {1 + a_1 z^{-1} + a_2 z^{-2}}
\tag{1}
$$

Has difference equation:

$$
y[n] = b_0 x[n] + b_1 x[n-1] + b_2 x[n-2] - a_1 y[n-1] - a_2 y[n-2]
\tag{2}
$$

How would you compute the difference equation when $H(z)$ is expressed as poles and zeros without expanding brackets:

$$
H(z) = g \frac {(z - d_0)(z - d_1)} {(z - c_0)(z - c_1)}
\tag{3}
$$

I'm happy to assume real poles and zeros for the sake of simplicity. I know that when performing a cascade on a computer that the output of the first section is fed to the second and second order sections are used. I'm interested in the maths though.

What I'm ultimately trying to do is get the difference equation for this expression direct without expansion:

$$
H(z) = g \prod_{k=1}^{K} \frac {b_{k0} + b_{k1} z^{-1} + b_{k2} z^{-2}} {1 + a_{k1} z^{-1} + a_{k2} z^{-2}}
\tag{4}
$$


texts: Phase response of S11/S22 parameter for an RF attenuator?

texts:
I am new into the RF measurements. I was measuring the S11/S22 (reflection) parameters of an RF attenuator (with different attenuation levels, e.g., 10 dB) via a vector network analyzer from Rhode & Schwarz.

I measured the magnitude of S11/S22 of the attenuator, which was typically below -30 dB. But when I measured the phase response of S11, then my teacher said it is not required.
So, my question is, why do we not need phase response of S11/S22? I understand that we seek linear phase response in the useful frequency range of S12/S21 parameters (i.e., constant Group delay).

Additional question: Do we not need it at all for other components as well, e.g., Directional coupler?

If not, then why?
If yes, then why?


texts: How can I improvise matlab code from GLCM to CLCM(Color Level Co occurrence Matrix) for feature extraction in texture analysis?

texts:
Currently I am doing my project on SKIN DISEASE DIAGNOSIS USING TEXTURE ANALYSIS in IMAGE PROCESSING domain.In this project,the histopathological image of diseased skin image is given as input and the name of the disease,details about the particular disease related to that image should be displayed as output.For feature extraction of such image,GLCM is not useful(as the input is a histopathological image).So,CLCM(Color Level Co occurrence Matrix) works well for such input.Can anyone help me with the matlab code for CLCM?

I am using MATLAB R2013a for this project.Is this okay or Do i need to change the version for better results?


texts: How is the luminance vector calculated for these linear matrix transforms

texts:
In the article "Matrix Operations for Image Processing" by Paul Haeberli there is an example described using a matrix to alter the saturation of an image. In the example given the author first converts the image to greyscale using custom values rather than ITU-R Recommendation BT.601 luminance vector.

Could anyone shed some light as to how he calculates those numbers?

  Where rwgt is 0.3086, gwgt is 0.6094, and bwgt is 0.0820. This is the
  luminance vector. Notice here that we do not use the standard NTSC
  weights of 0.299, 0.587, and 0.114. The NTSC weights are only
  applicable to RGB colors in a gamma 2.2 color space. For linear RGB
  colors the values above are better.


texts: Implemenation of Length N Moving average filter using block covolution algorithm(overlap and save algorithm)

texts:
I want to implement a Length N moving average filter using block convolution algorithm(overlap save algorithm)
I am copying block convolution algorithm(overlap save algorithm) from Proakis, DSP using MATLAB 3rd Ed
function [y] = ovrlpsav(x,h,N)
% Overlap-Save method of block convolution
% ----------------------------------------
% [y] = ovrlpsav(x,h,N)
% y = output sequence
% x = input sequence
% h = impulse response
% N = block length
%
Lenx = length(x); M = length(h);
M1 = M-1; L = N-M1;
h = [h zeros(1,N-M)];
%
x = [zeros(1,M1), x, zeros(1,N-1)]; % preappend (M-1) zeros
K = floor((Lenx+M1-1)/(L));         % # of blocks
Y = zeros(K+1,N);
% convolution with succesive blocks
for k=0:K
    xk = x(k*L+1:k*L+N);
    Y(k+1,:) = cconv(xk,h,N);
end

Y = Y(:,M:N)';                      % discard the first (M-1) samples
y = (Y(:))';                        % assemble output

Now using above function i am trying to implement a Length N moving average filter.
My script(m-file) is below
clc
clear all
close all
n=1:5
signal=n*pi;
subplot(221)
stem(n,signal);
xlabel('n');ylabel('amplitude');title('original signal')
noise=rand(1,5);
subplot(222)
stem(n,noise);
xlabel('n');ylabel('amplitude');title('noise signal');title('noise ')
p=ovrlpsav(signal,noise,6),
subplot(223)
stem(p);xlabel('n');ylabel('amplitude');title('corrupted signal')
m=5;
xu=0;
for i=1:m

    x=signal+noise;
    xu=xu+x;
end
xu=xu/m;
subplot(224)
stem(n,xu);xlabel('n');ylabel('amplitude');title('filtered signal')

I am getting a filtered signal similar to input signal as shown by above code,but am i following proper approach?


texts: Method for the Localisation of Impulsive Noise in Speech Signal

texts:
I'm sorry if this question has already been answered, but I couldn't find anything.
I'm working on a program for noise detection and need to locate the position of impulsive noise. The following image shows a part of the spectrogram of a file which contains a beeping noise at about 2 kHz.

I already tried analyzing the samples via the derivative method mentioned in this  paper, where the average derivative over the whole file is compared to that of the samples. The problem here is that the speech is very loud and no fast increase in volume can be seen.

I also tried to calculate the RMS for each window of the FFT magnitudes to create a threshold for measurement, but that does not seem to work either. Spectral leakage seems to make the detection quite hard despite using a window function.

I hope someone here can recommend a method to locate impulsive noise in high frequencies that I may have overlooked.


texts: Confusion regarding a library of Digital Signal Controller

texts:
I'll be using a NXP DSC with 56800Ex core and my aim is to implement a Elliptic IIR Filter in it.

The User Reference Manual of the 'General Digital Filters Library' (which is available here) has functions to implement IIR Filters.

But on page 3-19 shown above, it says

  The filter coefficients are calculated using Butterworth approximation.

Does this mean that we can only use Butterworth Filter? Also, as per my knowledge, we have to give the coefficients and therefore other filters can also be implemented.
So what exactly does the function do?
Can I implement a Elliptic filter using this function or not? If not, then what should I do use Elliptic filter?


texts: Phase of the sine signal by finding FFT of it in matalb using angle

texts:
I am using sine signal with  phase=0 an frequency=100KHz.
Computed FFT of it.
Using angle function and multiplied the result with 180/pi;
But getting 90.36 at 100KHz frequency.
Can anybody explain why I am getting 90 degree instead of 0 degree?


texts:
I am working on a seismic signal to extract long period long duration earthquakes. The sampling frequency is $2004\ \rm{ Hz}$ and the the signal has $20755200\ \rm{samples}$.  The spectogram seems to be very discrete and it is not smooth.

I need to work in windows in a range of $300\ \rm{sec}$, so when I zoom in, it is even more like blocks.

What I need to a smooth spectogram like this:

Here is the MATLAB code:

y=y-mean(y); % y is the input signal
fc = 8; % Cut off frequency
fs = 2004; % Sampling rate
[b,a] = butter(6,fc/(fs/2)); % Butterworth filter of order 6
x = filter(b,a,y); % Will be the filtered signal
[b,a] = butter(6,2/(fs/2),'high'); % high-pass filter
xx = filter(b,a,x);% high pass filter
xx   = sgolayfilt(xx, 6, 9); %smooth the data
specgram(xx,fs);
ylim([0 0.02]);
caxis([30 150]);

I appreciate your help.


texts: How to calculate the magnitude of rounding error in N-bit signed fixed-point arithmetic?

texts:
In example take into consideration 11-bit <11,5> signed fixed-point arithmetic.

  000000.11111 - 0.96875

  000001.00000 - 1

so the error magnitude will be:

  (1 - 0.96875)/2 = 0.015625 ?


texts: Estimating the power consumption of an arbitrary, real-time DSP algorithm

texts:
Consider a digital signal processing algorithm that operates over a block of $L$ [samples], with a real-time constraint $\tau$ [s] (i.e. its throughput must be $L/\tau$ [samples/s]).

The theoretical complexity of the algorithm is known, both in terms of arithmetic operations (additions, multiplications and comparisons) and memory accesses.

Is there any rule of thumb to get a (rough) estimate of the power consumption of this algorithm for typical hardware implementations (e.g FPGA, ASIC, General Purpose Processor, etc.)?


texts: What is the relation between hilb and Hilbert transform in MATLAB?

texts: In MATLAB, there is a command called hilb(n) which it says gives a Hilbert matrix of order n.  Is this matrix the one used for Hilbert transformation?


texts: Calculating similarity between two signal plots

texts:
I have two plots, each having frequency as x-axis and Gain as y-axis. By taking one data set as a reference I have to calculate similarity between them.

The graphs have same values on the x-axis and have same range on the x axis

Can 2D Correlation or Co-variance can do a decent job ? or
should I opt for Fréchet distance or DWT as I have read in some other posts ?

The first plot is the reference plot.

Here are the plots:

Please help!


texts: Periodic signal checking when using $\Sigma$

texts:
Is the following signal periodic?

$$ \sum_{κ=-\infty}^{\infty}\left[\mathrm{rect}\left(\frac{t+2κ}{10}\right)\right]+\cos\left(\frac{π}{75}t\right)  $$

where rect is the rectangular signal


texts: Moving average filter for sound in MATLAB

texts:
y[n] = b0x[n] +b1x[n-1]

I have this MA filter in matlab as following

 load handel;
 x = y;

 b0 = 0.5;
 b1 = 0.5;

 N = length(x); % Length of input signal
 y = zeros(size(x)); % Allocate space for outut

 y(1) = b0*x(1); % First sample, assuming x(0) = 0
  for n=2:N % Remaining samples
    y(n) = b0*x(n) + b1*x(n-1);
  end;

 soundsc(y,Fs);

I was asked to listen to the signal and describe it , well I have run the load handel without filter in matlab than i have run the signal with this filter , it sounds like that with filter the signal sounds more demped and less noise compared without filter and when i changed the b1 to b1 =-0.5 it hears only the descant the base of sound is gone . Now I want to know the reason behind.

What is the effect of the filter on the input signal which pass through the system?
what is the effect of the filter on input signal which pass through the system if we change b0= 0.5 and b1 = -0.5?


texts:
I have position and acceleration measurements and I would like to apply a Kalman filter to estimate the velocity of the system.

I am not sure yet about how to procede, but I will check the already answered questions on this website (like Estimating velocity from known position and acceleration , Kalman filter with accelerometer with DC offset and Applying Kalman filter to a data set).

Once I have understood how to procede, I will implement it using Matlab. There I saw there are 2 different types of Kalman filter: discrete and continuous. What is the difference? I mean, what is the difference of working with one or the other?

EDIT to be more clear, I am referring to function kalman and kalmd

Since I work with a set of data, should I use the discrete one?


texts: SDR : Signal divided between the edges of the receiving window

texts:
I'm using a SDR (specifically the RTL-SDR) to analyse the spectrum and i noticed a phenomena that does not seem logical.

The signal I want to capture is centred on 868.3 MHz and the central frequency of the SDR is centred on 869.3 MHz with a bandwidth of 2 MSamples. What I observed is that the 2 halves of the signal are displayed on both sides of the receiving window instead of just the half between 868.3MHz and 868.3MHz - 0.125MHZ ( as my signal bandwidth equals 125 KHz ) .

What do we call this phenomena and what is the explanation behind it ?

The first image is when the central frequency of the SDR is 868.3MHz

The second image is when the central frequency ofthe SDR is 869.3MHz


texts: Compare two signals, find times when they are affected by the same processes

texts:
I have two (noisy) signals $S_1(t)$ and $S_2(t)$. They contain information on three distinct events occurring in a physical system. Event $A$ affects only $S_1(t)$, event $B$ affects both $S_1(t)$ and $S_2(t)$, and event $C$ affects only $S_2(t)$. I want to identify the sequence of events in the system by analyzing the signals.

Is there some generic way to identify if two strings of data contain correlated events and find the timing of said events?

More specifically, my signals take the form
\begin{equation}
\textrm{d}S_i (t + \textrm{d}t) = S_i^{\textrm{true}}(t)\textrm{d}t + \textrm{d}W(t),
\end{equation}
where $S_i^{\textrm{true}}$ is the noiseless value of the signal calculated from the state of the system at time $t$, $\textrm{d}W(t)$ is a Wiener element (uncorrelated Gaussian noise), and the statement is to be read as an update formula. The l.h.s. is what is available to me. In case of an event the true values change sign in the following manner:
\begin{align*}
A &\Rightarrow\ S_1^{\textrm{true}} \to -S_1^{\textrm{true}}, \\
B &\Rightarrow\ S_1^{\textrm{true}} \to -S_1^{\textrm{true}} \textrm{ and }  S_2^{\textrm{true}} \to -S_2^{\textrm{true}},\\
C &\Rightarrow\ S_2^{\textrm{true}} \to -S_2^{\textrm{true}}.
\end{align*}
The way I am trying to do things right now: I am using a running average to get rid off some of the noise, then I find the second derivatives of the smoothed signals. If the true value changes, there will be a peak in the second derivative. I then determine whether there is a peak in both signals or just one at the given time. However, this does not work very well, as the derivatives pick up the noise, too.

Here is a plot of my raw and running window averaged signals:

I have no experience with signal processing. Is there a canonical procedure? If the question is too broad, where can I find recipes to deal with issues as such?


texts: Why do convolution kernels such as Gaussian, Laplacian, LoG almost always seem to be expressed in integers?

texts:
I'm a total newb in search of some deeper understanding, but I'm not able to read the maths behind these on Wikipedia.

If I understand correctly, you get the new value for each pixel by multiplying each of its neighbours' values by the corresponding member of the convolution kernel/matrix. For Laplacian there seem to be two commonly used such kernels:

 0 -1  0
-1  4 -1
 0 -1  0

and:

-1 -1 -1
-1  8 -1
-1 -1 -1

But why do they use only integers? I believe modern devices have floating point as fast as faster integer math these days.

Is this just a leftover from the past when floating point was slow and now it's become a kind of tradition? Perhaps people just re-using the old info they find from those days without deeper understanding? Or is there an actual reason why integers are used?

In fact, won't the use of integers here result in artefacting? What am I missing?


texts: Morphological operations to isolate and preserve network graph structure?

texts:
Suppose there is an image with long, connected lines which form connected, closed cycles in an image.
Is there a standard morphological approach to isolate only the lines, and drop all other pixels?

I can imagine a couple of slow approaches:

convolutional approach where a window -- say, 50x50 -- is passed along the image, and any non-zero pixels which are not connected to a pixel on the perimeter of the window are set to 0.

Skeletonize - threshold the components - close - skeletonize ...: the skeleton is calculated, small components are dropped, disconnected skeletons are closed, and the skeleton is recalculated and re-thresholded. Repeat.

But it seems like this would be a common problem, and there should be an easy solution...?


texts: How to implement Fourier Descriptor of an image?

texts:
I want to implement Fourier Descriptor of an object. I have read link. However, I have some questions about normalizing Fourier Descriptor.

First, if I want to normalize the position of the starting point, according to the above link, This is done by subtracting the phase of the second Fourier descriptor $\phi_1$ from the phase of all Fourier descriptors and weighting by $k$; that is:

$a(k)=a(k)e^{-j \phi_1 k}$

Is the weight $k$ equal to the $k$ in $a(k)$ ?

Second, how could I do in order to achieve rotational invariance when I implement Fourier descriptor ? Thank you very much for your reply.


texts: Is there a version of Welch's method that doesn't look for power?

texts:
Welch's method splits a time signal, $x(n)$ into $M$ periodograms $P_m$,

$P_{x_m,M }(k) = \frac{1}{M}|F_k(x_m)|^2$

and averages them to give the Power Spectral Density (PSD),

$S_{x}(k) = \frac{1}{K}\sum_{m=0}^{K-1} P_{x_m,M }(k)$

But if I am not interested in the power can I use something like an averaged Fourier transform to represent the Fourier transform for the whole signal $x(n)$?

I.e. The exact same as above but instead of averaging the periodograms, $|F_k(x_m)|^2$, I just average the fourier transforms $F_k(x_m)$ within each window?

Kind of like a STFT except averaging the windows, and is there any merit to using something like this?


texts: Difference Between `fft(signal, nfft)` and `fft(signal)`

texts:
I want to know the difference between these following commands (matlab) and why I get two different outputs?

fft(signal);

fft(signal, nfft);

I found that the output arrays are not the same, but I don't know why!
how it works with the signal of any length without padding the original signal with zeros!

signal length is't a power of 2
nfft is the next power of 2


texts: How to filter this signal $x= \sin(2\pi 10 t)+\sin(2 \pi 40 t)+\sin(2 \pi 50 t)+\sin(2 \pi 60 t)$ using comb filters in time domain?

texts:
Can anyone suggest a matlab implementation of this problem? I've done this uptill now,

fs = 4200;
t = 0.001:1/fs:1;
x= sin(2*pi*10*t)+sin(2*pi*40*t)+sin(2*pi*50*t)+sin(2*pi*60*t);
% figure,plot(t,x);
t1 = t - 0.05;
x1 = sin(2*pi*10*t1)+sin(2*pi*40*t1)+sin(2*pi*50*t1)+sin(2*pi*60*t1);
% figure,plot(t,x1);
yref = x - x1;
% plot(t1,yref);
N = 4096;
f = (0:N-1).*fs/N;
y_f = fft(yref,N);
plot(f,abs(y_f)),xlim([0 100]);
% freqz(y_f,1,1024,fs);

I am getting 10 and 50 Hz but from here how to remove 50 Hz?


texts:
I have just read on page 3 of this paper:

https://www.onsemi.com/site/pdf/Filtering_Hearing_Ai.pdf

That an "odd stacking" DFT can be created, essentially sacrificing the DC and Nyquist bands to offset the band centre frequencies by $Fs/(2N)$.

For the synthesis it describes the math as:

$ Y[k] = \sum^{n=0}_ {N-1} y(n) W^{-(k+v)n}, k=0,1,...,N-1  $

where v is 0 for an even-stacking DFT and 0.5 for an odd-stacking DFT.

So it appears to just be a $\pi$ phase offset on the complex exponential used on the DFT.

Have I interpreted this correctly? Can anyone explain further how this causes the offsetting of the frequency bands?


texts: Finding the impulse response, given a differential function involving the input and output

texts: I'm starting to teach myself some signal processing, and I'm stuck on this problem where I have to find the impulse response $h(t)$ of a causal LTI system given: $$y''(t)+2y'(t)+y(t)=x(t)$$ So far what I've done is taken the Laplace transform of both sides of this equation to get $$(s^2+2s+1)Y(s)-3sy(0)-y'(0) = X(s)$$ So the problem now is that I have all these inital values that I don't know what to do with. If I understand correctly, for all $ t<0$ in a causal system we have the output as zero, but this doesn't include zero. If anybody can point me in the right direction, I would really appreciate it. Thanks.


texts: How to determine STFT parameters as a function of desired resolution

texts:
None of the Q&As I've seen thus far have given a clear answer to what seems like a simple question: assuming I know the CenterFrequency, FrequencySpan/Bandwidth, SampleRate, and one or both of DesiredTimeResolution and DesiredFrequencyResolution, how do I determine the WindowSize and Number of Points for an STFT(nfft)? Working with complex/baseband signals if that makes a difference.

This answer seems to be saying that frequency resolution/resolution bandwidth is just SampleRate/nfft.
However that doesn't give me a formula for relating desired time resolution to STFT input arguments. Nor does it give me a formula for determining WindowSize (I've seen several explanatinos which say WindowSize should equal the NumberFFTPoints) but would like to double check that this is correct and/or understand what happens when they are different.


texts: Gesture recognition using only simple techniques

texts:
I'm following a computer vision course and I have this exercise: write a program that, given an hand image, it can recognize if the hand is open, closed, in a punch, or holding an "ok" posture, using only the techniques provided till now (pixel 4/8 connected, connected region, contour finding, holes finding, blob property like centroid, area, perimeter, eccentricity, image moments, image transformation like invert/power/log/gamma correction/contrast stretching, histogram computation and equalization).

I have done it with some basic blob properties(closed hand has a low eccentricity, "ok" has a hole, open hand has a big difference between the area of the inscribed ellipse in the blob and the blob area itself with a low eccentricity).. It seems to work but the first image is a little bit problematic.

I think there could be something more to make a more robust algorithm. Maybe some kind of moment property? Could some blob axis/orientation/extreme points help?

PS test images:


texts:
While searching for something FM, I came across this code. I am not able to understand a particular line of code. I have written a comment against that line of code.

fs = 1000;
fc = 200;
t = (0:1/fs:0.2)';
x = sin(2*pi*30*t)+2*sin(2*pi*60*t);
fDev = 50;
int_x = cumsum(x)/fs;  %%% This line of code
xfm = cos(2*pi*fc*t).*cos(2*pi*fDev*int_x)-sin(2*pi*fc*t).*sin(2*pi*fDev*int_x)   ;

I understand the cumsum function but I am not able to understand why is he dividing the integration by Fs i.e. sampling frequency.

Link to the code - https://in.mathworks.com/matlabcentral/answers/318969-fm-modulation-without-using-fmmod


texts: Smoothing "font characters" outlines in bitmap image

texts:
Adobe Acrobat has this as a feature called ClearScan. While doing OCR it allows image enhancement, by smoothing irregularities on font outlines.
For example, this source:

turns to this:

Are there other tools with such feature, or public algorithms (Matlab/Python) that can smooth contours of arbitrary font-like objects?


texts:
i have a sinus value at 80MHZ amplitude 0.5 with DC at 0.5 which was converted to digital sinus of values from 0-255.
I have peformed FFT and found exactly a DC tone at 127 and 127 tone at 80MHz.
both of them are pure digital values(not voltage)
SNDR is the ratio between the signal power and noise+distortion power.
how do i perform it in this case?

Thanks


texts:
i've read that it is possible to interpolate between two sets of filter coefficients (if they are close of course), but how to interpolate between two sets of coefficients which are of different orders? Is it possible ? If not what is the possible way to fade between two IIR filters ?

thanks

Jeff


texts: Steganography using wavelet transform — not every wavelet recovering the hidden message efficiently

texts:
I'm working on a project to hide simple string inside a 255x255 image using Wavelet transform. The problem is that,

When i used Haar Wavelet the result is successful,
This is what I've done,

Decompose the image into its approximation and details(Vertical,Horizontal,Diognal)
then Hiding message in the details of the Image,
Reconstructing the Image with Hidden message,
Applying same wavelet on the Reconstructed Image.
Decoding the algorithms applied for the encoding of the Message.

but when i tried these Wavelets,

Coiflets 'coif1', ... , 'coif5'

Symlets'sym2', ... , 'sym8', ... ,'sym45'

Discrete Meyer'dmey'

Biorthogonal'bior1.1', 'bior1.3', 'bior1.5'
'bior2.2', 'bior2.4', 'bior2.6', 'bior2.8'
'bior3.1', 'bior3.3', 'bior3.5', 'bior3.7'
'bior3.9', 'bior4.4', 'bior5.5', 'bior6.8'

Reverse Biorthogonal'rbio1.1', 'rbio1.3', 'rbio1.5'
'rbio2.2', 'rbio2.4', 'rbio2.6', 'rbio2.8'
'rbio3.1', 'rbio3.3', 'rbio3.5', 'rbio3.7'
'rbio3.9', 'rbio4.4', 'rbio5.5', 'rbio6.8'

I am getting weird characters rather than my original Hidden Message.
What does that mean ? How would i know which Wavelet should i suppose to use when the message is totally random ?


texts: Relaxation time in a stationary stochastic process

texts:
Foreword: I'll give some background here because I'm not sure all my premises are correct. I encourage you to IGNORE THE BACKGROUND, unless you're really interested / bored :)

The Question

I want to calculate the eigenvalues of this system.

$\frac{dM_0}{dt}=k+t_0 M_1 - M_0 (t_1+d)$

$\frac{dM_1}{dt}=t_1 M_0 - M_1 (t_0+d)$

In matrix form, I want the system to look like $m^{'}=Ax$, and I'd compute the eigenvalues from the square matrix A. But here I have $m^{'}=Ax+b$ -an inhomogeneous set of ODEs.

My intuition is that eigenvalues don't depend on $k$ ($b$) at all, but I have been unable to confirm this by any means (consulted many ODE/eigenanalysis texts, and the internet). So the question is: can I ignore $k$ when computing the eigenvalues?

Background

In order to construct the steady state probability distribution function (PDF), I run a simulation many times, record the "final" state of the system, and repeat. At the end of N simulations, I have a vector of final states $[X_1 X_2 ... X_N]$, which as $N\rightarrow \infty$ the histogram of which converges to the probability distribution function.

It seems to me that the optimal time $T^*$ to simulate each individual run is when the stochastic system has "relaxed" enough that it no longer "remembers" its starting condition. A good metric for this is autocorrelation: the lag time for which autocorrelation goes to zero (or some arbitrarily small value) seems like the time I should choose for $T^*$. If I choose a time $T<T^*$, every simulation is correlated with the start condition, and not independent observations. Thus, the PDF will never converge on the correct distribution. If I choose a time $T>T^*$ then I do "extra" computation. Since it turns out the simulation is VERY inefficient, even simulating a little past $T^*$ results in significant computational overhead.

The question is, how to choose $T^*$? Suffice to say, numerically computing it is out of the question. One idea suggested to me is to calculate the eigenvalues of the system. The inverse of an eigenvalue corresponds to the relaxation time of a system, which is equivalent to the lag-time tau at which autocorrelation = $\frac{1}{e}$. Thus, the time $T^*$ could be chosen instantaneously by selecting the longest relaxation time of the system (smallest eigenvalue)--or some multiple of this time to get the correlation sufficiently small. Bonus points if you weigh in on whether this is a good idea or bad, and why--but not necessary to answer the question.


texts:
Is it possible to represent an aperiodic signal using an array of N samples?
I am confused about this because you obviously have to window a function in the time domain to sample it.

Now what happens with aperiodic signals?


texts: What does "real time" signal processing mean?

texts:
I thought this was supposed to be an obvious question, until I finally set up my real time system.

So basically I have a transmitter that sends 128 samples/second to a receiver. The transmitted information is stored as an object in MATLAB and continuously updated.

When people talk about real time signal processing, I'm really confused as to what they mean by real time.

For example, say I want to extract the "mean" feature of this signal. Do I compute the mean when I receive one sample, two samples, all 128 samples, or ...what is this mean value?

More intriguing for me is the prospect of doing real time wavelet transform for joint time frequency analysis. Again, the question of "real time" comes up. How many samples do I need to reliable compute the wavelet (or fourier) coefficients so I could get a good view of the energy contained in this signal.

Can anyone who is knowledgeable on this topic please elaborate how when and under what condition do you compute "features" or perform frequency domain analysis for a real time system.

Thanks!


texts: Explain if the following processes are stable or not

texts:
I have got a question about my homework and was wondering if anybody could help. In the picture below you can see the processes. I need to determine whether they are stable or not.

Now I have only been able to get to this point (see next picture):

When a>0 the system should be stable. So looking back at the second process, this one should be stable too (because a = 10). However, when checking the anwsers, it says both the second and third process are not stable.

Could anybody explain how I can continue the second- & begin the third process?

With regards


texts:
I am wanting to compute the Z-transform of $f(n) = a^{|n|}$ . 'a' is a positive constant.

Looking at the transform table, I found that Z-transform for $a^n u(n)$
is available from the tables and is $\frac{Z}{Z-a}$. Where $u(n)$ is the unit step function.

I am trying to decompose $f(n)$ as --
$f(n) = a^n u(n)  +  a^{-n} u(-n) + \delta(n)$

Then using the table to find $Z( f(n) )$ as --
$\frac{Z}{Z-a} + \frac{Z^{-1}}{Z^{-1} - a} + 1$

Can anyone tell me if this approach is correct and/or suggest an alternate way.
Thanks!


texts:
I edited the initial question to make it general and thus applicable in more situations.

If the Fourier Transform of $x(t^2)$ is some function of $\omega^2$. Then what is the Fourier Transform of $x(bt^2)$, for $b>0$?

Is it $\dfrac{1}{b}$ times a function of $\left(\dfrac{\omega}{b}\right)^2$, or $\left(\dfrac{1}{\sqrt{b}}\right)$ times a function of $\left( \dfrac{\omega}{\sqrt{b}} \right)^2$, or neither?

Thank you for your time.


texts: How to generate sinusoidal images that shift over time

texts:
I try to generate an image using this code.

w=10*pi/256;
E=zeros(256,256,100);
for t=1:100
for y=1:size(E,1)
    for x=1:size(E,2)
        E(x,y,t)=255*(sin(w)*sin(w*y*2)+1)/2)+1;
    end
end
imshow(E(:,:,t),[])
title(t);
drawnow
pause(0.1)
end

However I'm confused how to generate the image that shift over time, like for time t=2, the stripes is shifted by 0.0001 in the x direction.

Any idea how to fix this? Thank you.


texts:
I want to use maximumu correlated kurtosis method described by McDonald to the vibration signal to extract bearing faults frequencies, but confused! what are the appropiate input parameters used in the matlab function, or how can I fix them:
[y_final f_final ckIter] = mckd(x,filterSize,termIter,T,M,plotMode)
where I used parameters that are alredy used on, while they gave me better results
filterSize=400

termIter=30

T=100

M=7

plotMode=1


texts:
What is meaning of multi band response?
I have already tried reading filters texts but no help.
How it is different from simple band pass or band stop


texts: What happens if I zero out the initial terms of the fft output and take inverse fft?

texts: I am new to signal processing domain. When I run fft on data of length 1000, I get 1000 complex numbers. Now, if I want to extract the low frequency information or signal approximation, I take the first few elements(say 50) of the array and also the same number of elements from the last(as it is symmetric). Excluding these values, i put everything as 0 and take inverse fft and I get back the approximation of the original signal. Now, having done the above, I take back the original array of 1000 complex numbers,  put the first 50 and last 50 numbers as 0, and then repeat the process on the remaining values to get successive higher frequency information. My main goal is that I want to separately extract low frequency and high frequency information from the signal(like in wavelet). Is my approach correct?


texts:
Let $\alpha, \beta > 0$ and $\Delta := \nabla^T \nabla$ be the discrete laplacian operator,
$$\nabla: \mathbb{R}^{n_x \times n_y \times  n_z} \rightarrow \mathbb{R}^{3 \times n_x \times n_y \times n_z}, \hspace{.5em}w \mapsto (\nabla_xw, \nabla_yw, \nabla_zw)
$$
being the discrete spatial gradient operator in $3$D. You can see $\Delta$ as a linear operator of shape $p \times p$, where $p := n_xn_yn_z$.

My intuition is that the operator $(\alpha I + \beta \Delta)^{-1}$ corresponds to some kind of Gaussian smoothing. Is this true ? If yes, what's the width of the corresponding kernel in terms of $\alpha$ and $\beta$ ?


texts:
So I was working on a program to comparing the similarity between the pronunciation of a word, First I got the audio file, and I use FFT on it to get the amplitude over frequency array, then I just square everything to get the Energy Spectral Density of it. The graph that can be drawn from those frequency of 4 sample of the same word I spoke myself look very similar to each other as shown in the picture I included here. All 4 sample have 2000 value, represent frequency from 0 to ~800Hz.

But since even when it's me saying the word 4 times, the frequency of the word cant be the same everytime I say that word, the graphs look visually similar but they are a few sample different from each other, like the 2nd highest point in the first graph is at the value number 1150, the 2nd graph have it at 1124, the 3rd have it at 1132 and the 4th have it at 1160, same thing happen with the first peak.

So my question is: how do I actually compare these 4 graph together to show how similar they are? Since what I was suggested to do the first time around was to using Mean Square Error but they didnt work well due to the slight different everywhere.


texts: How do I detect the sharp edges in the following image?

texts: I am looking to count the sharp pointy triangular things in the image. However I was unable to do so. One approach I did try was take a convex hull of the contour and then find convex defects but it led to lot of irregularities. You can find the it in the labelled image below.


texts: Is there any way to measure of Gaussian-ness?

texts:
I have some sampled data that has $1/f$ noise in it, with departures from the mean. These are long term departures. I could use something like a median filter but the window length would be longer than I would like.

Is there any way to measure its departure from a Gaussian distribution?
Or measure how 'Gaussian-like' a time series statstical sample is?


texts: Why can Fourier series be used on a non-repetitive function?

texts:
I was just wondering why Fourier series can be used on the function in the linked image. This is since I thought the function had to repeat itself to use Fourier series on it. Or is it saying a period is from 0-5s and so it will repeat every 0-5s? like the 2nd graph i made?

I can find the answer easily just not sure why fourier series can be used?


texts: Phase Information at Higher Frequencies in Continuous Wavelet Transform

texts:
I'm using the code I found here to compute the wavelet transform of a sine wave with a constant frequency.

#!/usr/bin/python2

from pylab import *
import matplotlib.pyplot as plt
import numpy as np
import scipy

x = np.linspace(0, 10, 65536)
y = np.sin(2 * pi * 60 * x)
N = len(y)
Y = np.fft.fft(y)

J = 128
scales = np.asarray([2 ** (i * 0.1) for i in range(J)])

morletft = np.zeros((J, N))
for i in range(J):
    morletft[i][:N/2] = sqrt(2 * pi * scales[i]) * exp(-(scales[i] * 2 * pi * scipy.array(range(N/2)) / N - 2) ** 2 / 2.0)

U = empty((J, N), dtype=complex128)
for i in range(J):
    U[i] = np.fft.ifft(Y * morletft[i])

plt.imshow(abs(U[:,scipy.arange(0,N,1)]), interpolation='none', aspect='auto')
plt.title("Sine Wave")
plt.xlabel("Translation")
plt.ylabel("Scale")
plt.show()

The result looks alright.

What I'm really interested in looking at is the phase. I'm extracting phase from the above code using:

imshow(np.unwrap(np.angle(U)), aspect='auto')

and it looks like this:

Why is the phase information present at frequencies higher than (or scales lower than, NOTE: inverted y-axis) that of the signal?


texts: Relation between histogram equalization and 'Auto Levels'

texts:
I need to enhance an image using only global linear adjustment of intensity, i.e.:

$$I'=aI+b,\qquad a,b\in \mathbb{R}$$

where $I$ and $I'$ are input and output (enhanced) image.

I already know how to determine best values of $a,b$ through histogram stretching given the minimum and maximum intensity within image. However, this will not work on the following image. I have added black and white pixel (in red circle):

The image already contains values 0 and MAX, further histogram stretching is not possible.

However, I tried "Auto Levels" feature in Paint.NET software  with the following result:

Looking on the dialog window, there is some extra factor that affects overall histogram shape (number in the middle right part):

Is this an additional scaling factor? Or it is some parameter of non-linear transform?

How is is determined?


texts: Pre-emphasis / De-emphasis filter design

texts:
I've designed a simple FM communication system and found that high-frequency noise is a real problem.  It appears that this is a well known issue, which has been solved through pre-emphasis of high frequencies at the transmitter and de-emphasis at the receiver.

The standard pre-emphasis filter appears to have the following form:

$$H_{pe}(f) = 1 + j2\pi f\tau_x, \quad \text{where } \tau = \rm75E-6$$

The de-emphasis filter is simply the inverse of that, or

$$H_{de}(f) = \frac{1}{1 + j 2\pi f\tau_x}.$$

I'm afraid I'm not much of an IIR filter designer, but from my research it appears to me that the de-emphasis filter is a relatively simple integrator, and that I can use the Impulse-Invariant method. The pre-emphasis filter, on the other hand, looks to be a bit tougher.  Should I use a bilinear transform approach?  When I tried to use the bilinear command in MATLAB it complained, presumably because there were no poles.

Any tips or hints would be greatly appreciated.

EDIT:
The MATLAB command that I used was:

[z, p, k] = bilinear([j*2*pi*75e-6], [], 1, 20e3);

MATLAB (using version R2011b) returned:

Error using bilinear (line 100)
First two arguments must have the same orientation.


texts: Are high pass/low pass filters FIR filters?

texts: I am trying to understand the basic categories of DSP filters, and it seems that the two main groups are FIR and IIR, however I'm not sure how to categorize high pass, low pass, band pass, and band stop (the 4 basic filters I learned about in university) and I can't find material on this either even though it seems so basic. Any help with making categories out of or distinctions with these groups would be greatly appreciated!


texts:
I would like to use a .WAV file and convert it into a spectogram. Any idea how to do this with Naudio? Basically the flow should be like this, I think:

Convert WAV file to PCM stream
Convert PCM stream to double-array and applying FFT
Using the result from (2) and convert it to a spectrogram


texts: The process of my DSP class used the Oppenheimer DSP book as text. Although I can easily understand the math in the book, I found it lack examples and applications, hence I still cannot build a picture how people actually use those under different circumstances. Can someone recommend me a DSP book that contains more practical examples and applications more than "This example uses the mathematical formula above to derive some numbers from some numbers"? Thanks in advance.


texts: How can it be that there is a series of integrals in Fourier Series if it’s a projection on a continuous basis?

texts:
If the process of finding Fourier coefficient is finding the projection of a signal on a member from an orthonormal basis, basis which is continuous in frequency. How can it be that Fourier coefficient is an integral and signal representation, which is the sum of these coefficients multiplied by the projections, is a series? The coefficient basis members are continuous, how can they be discrete in the signal representation while being continuous in the coefficient representation? The logic says: we took projection in continuous space, then we shall continue in continuous space while we multiply by that projection. Why do we have a mix of integral and a series then?

  How can it be that there is a “series of integrals” in Fourier Series if Fourier is just an operation of projecting a continuous basis and multiplying the projection on the (continuous) basis members once again?


texts: How to adapt 4th order IIR filter coefficients from matlab to cmsis dsp biquad struct

texts:
I'm designing an IIR 4th order filter to be implemented using CMSIS DSP on a cortex-M4 using a direct form II trannsposed structure. I get this coefficient from matlab's "butter" function:(consider Fs = 48kHz, bandpass Fc1 = 445 Hz, Fc2 = 561Hz)

b[5.70278412657461e-05   0   -0.000114055682531492   0   5.70278412657461e-05]
a[1   -3.97002047467733   5.9189048757951   -3.92762130007361   0.97875498919531]

I know i have to use two stages of a biquad cascade IIR filter and the according structure should receive coefficients in the follow way:

float32_t[10] = {b10, b11, b12, a11, a12, b20, b21, b22, a21, a22};

Question: Can some tell me how to translate matlab answer into CMSIS struct format?.


texts: Frequency of the wave in frequency domain

texts: If we have a 1-dimensional wave in time domain, it can be represented in frequency domain with x axis indicating the frequency of the wave and y axis indicating amplitude/magnitude of the wave. But what about a 2-dimensional wave in case of images in Frequency Domain? Because in 2-dimensional wave we have a cordinate (x,y) to denote the magnitude. So what is the frequency in this case?


texts:
I have signal which contains some noise only regions (glitch with the recording it seems).

The question is:

I would like to identify the discrete time samples where this noise part starts and ends. There are about 3 to 4 seconds of this noisy part available, there are 4 to 5 such instances, the total recording length is 153 seconds. The sampling frequency is 48 kHz.
I would also like to get to characterize the spectra of noise so that I can possibly try and remove it. How can I do this?


texts: What is the relationship between scales and frequency in a Morlet wavelet transform?

texts:
I'm using PyWavelets, with a complex Morlet wavelet. Its complex wavelet transform function requires scales as one of its parameters, rather than frequencies. However, I don't understand the relationship between the two. I'm generating my scales like this:

base_scale = CENTER_FREQ * quality / base_freq
far_scale = base_scale / 2**num_octaves
scales = numpy.geomspace(base_scale, far_scale, num=num_octaves*voices_per_octave+1, endpoint=True)

That is, if I expect my signal to be at about 20 kHz, I want to set the scales to be such that the frequencies returned are between 10 kHz and 40 kHz. In this case, the num_octaves variable might be 2, base_freq would be 10000, and voices_per_octave might be 50. Then scales would be an array of length 101, ranging logarithmically from the scale corresponding to 10 kHz up to the scale corresponding to 40 kHz.

However, this only works if CENTER_FREQ (the wavelet's central frequency) is a small number. If I increase it beyond 31, the frequencies returned by:

coeffs, freqs = pywt.cwt(data, scales, wavelet, 1/quality)

(where quality is the sampling rate of the data) are dead wrong. Is there any easier way to determine the needed scales that correspond to the frequencies I want?


texts: Will a Bluetooth 5.0 USB dongle implicitly support aptX?

texts:
Am I asking this in the right place? I also asked on Super User and haven't got a response yet, but here's my issue...
I have an aptX low latency receiver/transmitter 2-in-1 which uses a 3.5mm audio jack, so I have that transmitting audio from the headphone jack of my device. Now, I want that audio to show up on my Windows 10 PC. I could just get another of those modules in receiver mode and plug into the microphone jack on my PC, but that compromises sound quality and for some reason there's a massive 100ms delay on it which defeats the purpose of aptX-LL so I need a dedicated USB receiver.
I can find every aptX USB transmitter under the sun, which will be useful for a different part of the project, and I can find lots of USB receivers that don't mention supporting aptX, but I can't for the life of me find a receiver that I can plug into my USB port, and have it show up as an audio device and get that critical sub-40ms latency... At least not for less than a $70 command station, which seems excessive considering everything else including the transmitters seem to be in the 30 range so... What am I missing?? If a regular USB receiver is Bluetooth 5.0 does that mean it implicitly supports aptX? I'm so confused why this doesn't seem to be a thing that exists!


texts: What is the minimum number of time points recommended for accurate instantaneous phase estimates using the Hilbert transform

texts: I have two time series.  The two series look like they follow the same function, but one appears to lag slightly behind the other.  I would like to determine if there are phase differences between the points in these time series, but each time series is only 4 points long.  Is this enough data to get accurate phase differences using the Hilbert transform?


texts: What is correlation kernel and compare with gaussian kernel

texts:
I read a paper that said about correlation kernel that defined:
$$W(x-y)=(α/1+d(|y − x|))$$
where $α =  (\int(1+d(y − x)dy)^{-1}$, $(d(|y − x|))$  is spatial Euclidean distance from the central pixel.
My question is what is correlation kernel and benifit of that kernel comparison with gaussian kernel? Thank you so much.


texts: Terminologies - sparse channel, sparse input, compressed sensing

texts:
The term sparse in general means that there are more elements that are zero valued or very close to zero in comparison to the number of non-zero.  In speech deonvolution research papers, the channel is assumed to be sparse, so  the channel has more number of zeros. For instance: the channel is sparse when the speech signal is recorded by one microphone in an enclosed space.  The microphone signal is corrupted by reverberation and additive noise.  I don't understand why the channel is sparse in this case.

In compressive sensing methods, I have seen that the input signal is sparse. In this case, it is related to nyquist sampling theorum and there should be a large number of measurements. I have the following questions:

Under what conditions is the channel sparse and what is the significance of sparse?
Under what conditions is the input information sparse and why?
When is the term compressed sensing used? Is it applicable to sparse channel or sparse input?

Please correct me where wrong.

I could not find definite answers and information in research articles and book related to these questions, and shall really appreciate a succint explanation. Thank you


texts: Modelling room impulse by raytracing method

texts:
I'm building a raytracing engine to get an approximation of the impulse response of a room, but the problem is that I don't really understand what to do after the ray tracing.

After I get all the reflections that collide with the receiver, I get all the rays with their acoustic energy and collisions time. But how do I turn this into an actual impulse response that i can hear using convolution. I was thinking in using a three band eq and turn the reflection energy into intensity, to set parameters for the eq, then send an impulse through the filter and do this for each ray; and finally sum all the impulse responses at their correct times. Please help, I'm not sure if the method I'm going to use is correct.


texts: How to detect decaying oscillations in a signal

texts:
I have a sampled signal which contains (among other things) some damped oscillations. From the underlying theory I have reasons to believe that it can be characterized by some specific frequencies of oscillation $\omega$ and rates of decay $-\sigma$. I believe the signal contains components of the form $f(t)=e^{\sigma t+i \omega t}$. It might be that these are "there" only in some part of the signal, because at late times other components win over (due to exponential damping of these oscillations). I am interested in finding these dominating timescales. I need a hint how to proceed. There are plenty of tools available, like Fourier analysis, wavelets etc. However, I don't know which is the most appropriate for this task, as I have no experience with any of them, except maybe the Fourier transform. Below I include a graph of the signal that I want to analyze.

http://i.imgur.com/fZwhdkq.png


texts: Output of MATLAB's two-dimensional DWT function $\tt dwt2$

texts:
This is a conceptual doubt about the two-dimensional DWT. I am trying to understand the decomposition step of a two-dimensional DWT.

In the MATLAB explanation, the two-dimensional DWT generates 4 new images (cA, cH, cV, and cD) that represent the approximation, the horizontal details, the vertical details and the diagonal details.

I am also reading the explanation about this two-dimensional DWT. But the author shows the results as $\text{LL, LH, HL and HH}$.

I want to know if it corresponds to the same thing as cA, cH, cV, and cD. Is the following comparison correct?
\begin{align}
&\texttt{cA} =\text{LL} =\text{approximation}\\
&\texttt{cH} =\text{LH} =\text{horizontal}\\
&\texttt{cV} =\text{HL} =\text{vertical}\\
&\texttt{cD} =\text{HH} =\text{diagonal}
\end{align}

I'm also a little confused about the $\text{LH}$ and $\text{HL}$ results. In some posts, I have read that the $\text{LH}$ represents the horizontal details and in other, I have read that it represents the vertical details.

For example, as answered here , the user explains that:

  In the wavelet terminology, LL is the approxmation image, LH is the vertical details, HL is the horizontal details, and HH is the diagonal details.

But in this article the author says:

  The LL band corresponds roughly to a down-sampled (by a factor of two) version of the original image. The LH band tends to preserve localized horizontal features, while the HL band tends to preserve localized vertical features in the original image. Finally, the HH band tends to isolate localized high-frequency point features in the image.

Can someone clarify this for me?


texts:
In the literature, I see "corner detector" and "feature detector" are interchangeably used for these methods e.g. Harris, Shi-Tomasi, FAST, etc.
What is the difference between a corner detector and a feature detector? Are corner detectors (Harris, Shi-Tomasi, etc.) one kind of feature detector?


texts: SIFT between 1 image and dataset (2 images)

texts:
I use SIFT, and I want to find the angle of inclination of an image compared to two reference images, in short:

I have 3 images :

Image x: input image

Two images 1 and 2, reference images that are very similar.

My goal is to find the angle of rotation between the image x and a combination between the two images 1 and 2. My problem is that I dont know how to combine the two images of the database to compare them with to the image x.

The purpose of the combination is to extract a vector Keypoints and Mat descriptor from the 2 keypoint vectors and the 2 Mat descriptors.


texts: What are True Positive and False Positive with reference to image processing

texts: Say for instance, face detection algorithms - I see researchers testing their algorithms and table their results under TP and FP. What are these?


texts: Frequency contour detection in spectrogram

texts:
I have an image of a spectrogram and I wish to detect the tracks/contours of prominent frequencies present in the spectrogram.

In the end, I want to be able to get various prominent curves from the image.
The second image represents the expected frequencies. I want to be able to detect those frequency profiles are present in the noisy spectrogram. At the end, I expect a yes or a no as an output.

PS: Clearly I can see and tell that the expected frequency profile is present in the spectrogram. I wish to automate the process using a program preferably in python.


texts: what is the relationship between a spectrogram and the uncertainty principle heisenberg?

texts:
what is the relationship between these two things
Perhaps more resolution in a spectrogram is equivalent to knowing more the position of the electron and less resolution is knowing the velocity of the electron.


texts: Signal power and signal-to-noise ratio with repeaters

texts:
I am trying to solve this question but I am confused on some aspects

When solving part B what power signal will we use. Will we use power signal before repeater 2 or after repeater 2 (x5)
When finding maximum bit rate what should the SNR value? should it be average of all the SNR or what?


texts: Initial rest condition applied on $x(t)$ vs $h(t)$

texts:
Define the LTI system $\mathcal{H} : x\mapsto y$
Define the convolution for continuous-time system :
$$
(x*h)(t)=\int_{-\infty}^{\infty}x(\tau)h(t-\tau)\;\text{d}\tau
$$
The initial rest condition states that :

An LTI system with an input signal $x(t)$ is causal if and only if $x(t)=0$, $\forall t<0$

Now I have noticed that in some textbooks, the definition of convolution would have $h(\tau)$ instead of $x(\tau)$ and $x(t-\tau)$ instead of $h(t-\tau)$ (since $x*h=h*x$) but the initial rest condition now becomes :

An LTI system with an impulse response $h(t)$ is causal if and only if $h(t)=0$, $\forall t<0$

My question is :

Why does commutativity allow that for a LTI system to be causal then in one definition $x(t)=0$, $\forall t<0$ while in another definition $h(t)=0$, $\forall t<0$ knowing that $x$ and $h$ are two different terms?


texts: How is $-6$ dB pressure amplitude is the same as FWHM?

texts: I keep reading many literature, and in most cases $-6$ dB is referenced as the full width at half maximum (FWHM). Mathematically, how is that possible from the definition of dB?


texts: Synthesizing harmonic tones with Karplus Strong

texts:
I'm building my own MIDI synth.

I have to create a decent sounding tone for every note on the piano keyboard (A-1 through C7, or MIDI 21 through 21+88=109)

I've been using a basic ring resonator, where I create a ring of the right size (for example, sampling at 44.1 kHz the ring for A440 a.k.a. concert pitch would be 44100/440 i.e. 100.25 ~ 100 samples), fill it with static (values between -1 and +1) and then walk through it, averaging consecutive values.

The problem is that A-1 rumbles for 30 seconds, whereas C7 dies almost instantly.

I never noticed this problem before as I was working with a much smaller range; < two octaves.

Attempting to balance out the decays by multiplying by a suitable decay constant, say X[n] = .995*(X[n-1] + X[n-2])/2, where I choose the constant carefully for each note, is also failing.

Problem here is that even with decay=1.0, a high note still disappears almost instantly.  And it makes sense thinking about it; if the ring is only 10 samples,  even in 1/100sec we have covered 440 samples i.e. 44 revolutions, which will have smeared everything pretty close to 0.

One method I tried is feeding energy into the ring; creating a buffer of static enveloped by an exponential decay, and feeding this buffer into the ring.   It isn't working very well.

What I found works better is linear interpolation, so instead of doing ( X[n-1] + X[n-2] ) / 2, I do maybe:

0.98 * X[n-1]  +  0.02 * X[n-2] for a high note,  and
0.6 * X[n-1]  +  0.4 * X[n-2]  for a note lower down

etc.

But now the question is: how to choose the interpolation factor depending on the note?

I need a much better model to work with.

I think I also need to implement fractional delay;

Samp[ N ] += k *  Samp[ N - L ],

and if L = 10.2 say, I would have to do:

Samp[ N ] += k * ( 0.8 * Samp[10]  +  0.2 * Samp[ 11 ] )

so that is a completely different linear interpolation.

So my main question is: How to mathematically determine the necessary interpolation factor for a given note, so that all notes have an equal decay?

Supplementary questions are:
 - do I need to use a decay constant in addition to adjusting the interpolation?
 - is my understanding of fractional delay sufficient to get a working result?

If anyone knows of any resource for simulating a harmonic tone that is reasonably accessible, please do put it in as an answer. I've looked through JOS's resources, but he writes at a very advanced level, and I have difficulty understanding the material.

π

PS I can't find an appropriate tag;  pitch, note, gen*, tone, synth*, none of these are available.


texts:
I have a quadrature (IQ - 0..255) data byte stream that I receive from an rtl-sdr dongle as IQIQIQIQ.....

That means that negative values are < 128 and positive values >= 128

I am attempting to decode a 2FSK signal:

To normalize the data I do the following and get a resulting complex number with values ranging between -1 and 1 ... (more or less)

 Real = (I - 127.4) * 1/128
 Imag = (Q - 127.4) * 1/128

I demodulated with GNU Radio:

The scope plot correspond to the data that I get from a CC1101 tranceiver with an Arduino:
00000000011101110100010001110111010001000111011101000100010001110100010001000111011101000

This is the parameters for the CC1101 (SmartRF Studio)

The project is for a USB receiver to receive and identify signals from pager buttons. Unfortunately the CC1101 only works well for about 60% of the pagers.
That is why I want to use a RTL-SDR dongle and do my own demodulation.

GNURadio is great but it is a few hundered MB there is no minimal install ;-(

I would like to use the Goertzel algorithm to demodulate but it requires Real(float) input (not complex). If this is possible I could translate and compile the source to C or Delphi and create an executable that is maybe 1MB is size and a lot easier to deploy.

How do I prepare/use this data for the Goertzel algorithm?
Here is the prototype of the Goertzel C# code posted in another thread, that I want to use.

float goertzel(int numSamples,int TARGET_FREQUENCY,int SAMPLING_RATE, float* data)

It requires the data as an array of float. float* data
https://dsp.stackexchange.com/questions/23749/normalizing-the-magnitude-of-goertzel-filter-c

From reading the code I can see that it is not an array containing complex numbers (just an array of float) and I don't know how to get my IQ data into the format required by the Goertzel algorithm.

Finally my appologies and thanks to Marcus Müller (answered below) who spent a lot of time and effort to answer my original post in great detail, which had some incorrect, incomplete and confusing information, which I have now updated to be more complete.


texts: Implementing FIR/IIR "alike" filter in R

texts:
I am new to signal processing. I am trying to simulate something similar to IIR/FIR filter with $k$ delays to imitate acoustic echo reflection. The difference equations for FIR and IIR respectively are as follows:

\begin{equation}
y(n) = x(n) + \sum_{D=1}^kA(n)x(n-D)+v(n)\;\;\;\;\; (1)
\end{equation}
\begin{equation}
y(n) = x(n) + Ay(n-D(n))+v(n)\;\;\;\;\;\;\;\;\;(2)
\end{equation}
where $D$ is a delay in samples, the coefficient $A(n)$ describes
the changing attenuation related from object reflection and
$v(n) ∼ N(0, 10^{−3})$ is the noise.

  Equation $(1)$ and $(2)$ can be found on section $VII$, second
  paragraph of this paper

How to can I implement this? I started by writing the following code in R.

Edit after the answer

install.packages("signal")
t <- seq(0, 1, len = 100)
x <- rnorm(100) + rnorm(length(t),0,0.001)
y <- filter(Arma(b = 0.1, a = 0.1), x)

Unfortunately, the above approach does not allow me to have $a=0$ or $b=0$.

Remark

Maybe following function is generating equation $(1)$:

fir1(39, 0.3)+rnorm(40,0,0.001)

The second equation perhaps is called flange IIR filter, where
the delay is not constant, but changing with time. This effect
imitates time stretching of the audio signal caused by moving
and changing objects in the room.

Response to the answer below

The ARMA equation

$$y(n)=\sum_{0}^Ma(m)x(n−m)+\sum_{k=1}^Kb(k)y(n−k)$$

i.e. AR part is as follows:

$$y(n)=\sum_{0}^Ma(m)x(n−m)$$

and MA part is as follows:

$$\sum_{k=1}^Kb(k)y(n−k)$$

which don't match with equations $(1)$ and $(2)$. In equation $(1)$ the coefficient $A$ depends on $n$. The coefficient update at each iteration.

Probably an answer

Parameters

t <- seq(1, 4000, by = 1)
x<- sin(2*pi*t*2.3)
A<- rnorm(4000)
v<- rnorm(4000,0,0.001)
k <- 20

Equation $(1)$

for(i in 1:4000){
  if (i>k){
    y[i]<- x[i]+ A[i]*sum(x[(i-(k-1)):i]) + v[i]
  }else{
    y[i]<- x[i]+ A[i]*sum(x[1:i])
  }
}

Equation $(2)$

for(i in 1:4000){
  if (i>k){
    y[i]<- x[i]+ A[i]*(y[i-i%%k]) + v[i]
  }else{
    y[i]<- x[i]+ A[i]*(y[1]) + v[i]
    }
}

Does the solution make sense? Or do I need to generate $x(n)$ using FIR and IIR filter?


texts: Why the requirement of the GCD of the lengths of all circuits in the graph being one?

texts:
I am reading A Mathematical Theory of Communication. The second requirement of an ergodic process confuses me (emphasis mine):

  All the examples of artificial languages given above are ergodic. This
  property is related to the structure of the corresponding graph. If
  the graph has the following two properties the corresponding process
  will be ergodic:

  The graph does not consist of two isolated parts $A$ and $B$ such that it is impossible to go from junction points in part $A$ to
  junction points in part $B$ along lines of the graph in the direction
  of arrows and also impossible to go from junctions in part $B$ to
  junctions in part $A$.
  A closed series of lines in the graph with all arrows on the lines pointing in the same orientation will be called a “circuit.” The
  “length” of a circuit is the number of lines in it. Thus in Fig. 5
  series BEBES is a circuit of length 5. The second property required
  is that the greatest common divisor of the lengths of all circuits in
  the graph be one.

The author later elaborates

  If the first condition is satisfied but the second one violated by
  having the greatest common divisor equal to $d > 1$, the sequences
  have a certain type of periodic structure. The various sequences fall
  into $d$ different classes which are statistically the same apart from a
  shift of the origin (i.e., which letter in the sequence is called
  letter 1). By a shift of from $0$ up to $d - 1$ any sequence can be
  made statistically equivalent to any other. A simple example with $d =2$
  is the following: There are three possible letters $a, b, c$. Letter
  $a$ is followed with either $b$ or $c$ with probabilities
  $\frac{1}{3}$ and $\frac{2}{3}$ respectively. Either $b$ or $c$ is
  always followed by letter $a$. Thus a typical sequence is

  $$a\,b\,a\,c\,a\,c\,a\,c\,a\,b\,a\,c\,a\,b\,a\,b\,a\,c\,a\,c.$$

  This type of situation is not of much importance for our work.

But unfortunately, I still didn't get his point. Can you show me how to perform "a shift of from $0$ up to $d - 1$" to make two sequences statistically equivalent? In fact, the process already seems ergodic to me, but that's just my intuition, and I cannot justify my argument with rigorous reasoning.

In addition, what makes the above-mentioned process not of much importance in communication?


texts:
I have a question regarding an implementation of a modulator.
I guess some people here have an FPGA background.
I would like to implement a QPSK modulator of  31 Mbit bitrate using a 125 MHz clock.
I plan to generate the bitrate clock using MSB with a 32 bit phase accumulator also what can be called an NCO.
A 125 MHz clock has a 8 ns period.
A 31 MHz bit clock has a 32.26 ns period.
My bitrate clock will have an instantaneous period of 32 or 40 ns (but a mean of 32.26ns)
Is this a big issue to generate that kind of jitter on a modulator?
Is there a performance impact on a receiver?
Or maybe there is an another way implementing this?
Thanks


texts: I am sampling two sinusoidal analogue signals that are phase separated by $\frac{2\pi}{3}$, and there appears to be a significant amount of common mode noise in my samples. What, if any approach can I take to improve my SNR?


texts:
I'm trying to use and understand SciPy's deconvolve for a project I'm working on. I'm having some trouble understanding how to use it.

What I would like to do is to take two PMFs from discrete gaussian distributions and recover an unknown distribution using deconvolution. My understanding is that if I deconvolve the PMF from ~N(10, 1) and the PMF from ~N(30, 2), I should recover the PMF for the distribution ~N(20, 1). How can I do this in SciPy. It would be great, if I could see an example. Thanks in advance for the help.


texts: Hilbert FIR Filter for Q - Matching Filter for I

texts:
For an ultrasound application, I am collecting data from an ADC using high spec Altera FPGA which grabs data and can then do some DSP processing that would be otherwise too slow for a desktop to do.

Part of the DSP required is the precursor to envelope detection. The incoming signals from the ADC contain only the in-phase component of the signal. As a result I also need to generate the quadrature component to be able to then do the sum of squares of I/Q to get the envelope. To get the Q component a Hilbert FIR filter is called for.

Using MATLAB I've been able to design a suitable Hilbert filter for the application. However there are limited resources, so the filter is limited in terms of taps. Using the "Altera FIR Compiler II" tool I've then generated the filter for the FPGA.

I've been doing some simulation and the filter does correctly produce the quadrature component - for example if I feed in a sine I get a cosine. The trouble is the filter causes an attenuation of the signal - i.e. it isn't unity gain. Part of this will be due to coefficient rounding, and part of it due to the calculations being done in fixed point.

Using this filter if I try to do envelope detection using sum of squares, the resulting envelope is not brilliant. in the sine/cosine example I would expect the result to be a DC level, however it is in fact a DC level with a sine wave floating on it. This is presumably due to the Q signal being slightly attenuated. In fact if I multiply all the I values by a constant scalar found by trial and error, the envelope improves in accuracy.

The trouble is the pass band of the Hilbert filter isn't perfectly flat, there are <0.1dB ripples in it - there is a limitation on the number of taps I can implement due to resource restrictions. As a result the envelope of signals that are not pure tones get distorted again because my scalar for the I data is perfectly flat (its a simple multiplication).

What I am wondering then is can I improve the envelope by using a FIR filter on the I data as well. Rather than just delaying the I data and multiplying it by a constant scale factor, is it possible to design a FIR filter whose pass band response is matched to the Hilbert filter? If so, is there a simple way to do it using MATLAB?


texts: how to determine a length of a channel? ( given path, delay)

texts:
Given a channel matrices for each subcarrier of a frequency selective channel (OFDM): i.i.d..
channel =(randn(ntx,nrx, N)+1i*randn(ntx,nrx, N))/sqrt(2); simulated via channel simulator and a number of paths, delay

With these parameters, Can I define a length of the channel or a number of clusters $L$?


texts: What does 'correlation' and 'de-correlation' mean in image-processing?

texts:
Here are some problems:

Does it mean these is some relationship among neighboring pixels?
And I've read about DCT has a property of de-correlation,I have a doubt whether it means DCT neighbor coefficients have no relationship?
I've read some paper which says diffrence 2D-array(DCT coefficient subtracting by its neighbor) can weaken the relation of coefficients.How to understand this?

I would really appreciate if someone can help me!


texts: How to prove Line Spectral Pair coefficient properties?

texts:
I am trying to understand the definition of LSP coefficients from http://www.dspcsp.com/pdf/lsp.pdf.
We take the polynomial $a(x)$ of degree $M$ from the denominator of the LPC system function, and represent it as a sum of a palindromic and an anti-palindromic polynomial:

$a(x) = 0.5(p(x) + q(x))$,
where

$p(x) = a(x) + x^{M+1}a(x^{-1})$

$q(x) = a(x) - x^{M+1}a(x^{-1})$

Then, the textbook says,

  it is not hard to show that if all the zeros of a(x) are inside the
  unit circle, then all the zeros of $p(x)$ and of $q(x)$ are on the
  unit circle. Furthermore, the zeros of $p(x)$ and $q(x)$ are
  intertwined, i.e., between every two zeros of $p(x)$ there is a zero
  of $q(x)$ and vice versa.

How can this be proved?

UPD.

I found an article by Soong and Juang, "Line spectrum pair and speech data compression", with a proof. It depends on the equality
$ |z_iz - 1|^2 - |z-z_i|^2 = (1-|z|^2)(1-|z_i|^2) $,
where neither $z_i$ nor $z$ are assumed to be real.
But then the equality does not hold:
$ |z_iz - 1|^2 - |z-z_i|^2 = (z_iz - 1)(\overline{z_iz} - 1) - (z - z_i)(\overline{z} - \overline{z_i}) = |z_i|^2|z|^2 + 1 - zz_i -\overline{zz_i}  - |z|^2 - |z_i|^2 + \overline{z}z_i + z\overline{z_i} = (1-|z|^2)(1-|z_i|^2) - (z - \overline{z})(z_i - \overline{z_i})$


texts: Would a visible light beamforming "telescope" be feasible?

texts:
If you had a phased array of high speed photodiodes that were spaced sufficiently far apart, would it be possible to apply beamforming theory to create a steerable visible light "telescope"?

Traditional, earth-bound telescopes are limited by the size of their primary mirror due to the mechanical properties of figured glass as the size increases.  Additionally, the larger the primary mirror the more robust the tracking mechanism needs to be to handle the weight.

It seems that if you could instead "virtually" track the object by beamforming the output of a horizontally mounted photodiode array, you are no longer limited by the size of a primary mirror and/or infrastructure to physically move the telescope.  Angular resolution could be increased simply by making the array larger.  Furthermore, your tracking accuracy is only dependent on your software processing and electronics, not moving parts.

Of course, the problem likely becomes dealing with immense quantity of data collected from the photodiode array, and the sheer physical size of the array.

Light travels at about 299,700 km/s in air, so ~29.97 cm/µs.  That seems within the realm of high-speed electronics given a large enough array (10x10m?)

Difficult, huge amount of data...but it doesn't seem impossible unless I'm misunderstanding how beamforming works?


texts: Plotting Frequency Response Magnitude and Phase for first order all pass filter

texts:
Having trouble plotting the frequency response characteristics for first order all pass filter. The Magnitude is expected to be constant across entire freq and Phase is expected to be only decreasing (with jump from -pi to pi). I don't see either of this happening, provided plot from octave.
$H_{ap}(z) =  {z^{-1} - a^* \over 1 - az^{-1}}\tag{1}$
Above is condition to enable all pass filter.
code.m
num = [0.8,1]
den = [1;0.8]
TF = tf(num,den)
[zz,pp,kk] = tf2zp(TF)
figure(1);
bode(TF);
figure(2);
pzmap(TF);

Octave output:
Octave output:

Transfer function 'TF' from input 'u1' to output ...

      0.8 s + 1
 y1:  ---------
       s + 0.8

Continuous-time model.
zz = -1.2500
pp = -0.80000
kk =  0.80000


texts: Entropy : do we prefer higher or lower entropy?

texts:
The entropy of a source $H(S)$  gives the average codeword length to encode a given source alphabet. i.e. it is the average number of bits per symbol required to encode the information in the source. Let, $S_1 = 101101...$ be a message of length $N$. The entropy of the system that produces the message depends on the number of unique symbols, $k$ and the formula is $H(S) = ln k$ where $k$ is the number of unique symbols. The entropy $H(S_1) = \log(2) = 0.69$ bits/symbol. For another different message, $S_2$ having $k=8$ unique symbols, $H(S_2) = \log(8) = 2.07$ bits/symbol. I have 2 questions

Is there a way to determine $n$, the minimum word length for a message sequence of $N$? For example, $S_1 = (10)(11)(01)$ here $n$ = 2.
It is desired that the entropy of the system be maximized, so $H(S_1) < H(S_2)$ so should I prefer $k=8$ symbols then? Which entropy do we prefer - message with higher entropy and hence is this the reason why we prefer higher modulation?


texts: What is the link between equation of a continuous signal versus equation of its sampled form?

texts:
I have just started DSP self-learning. I am a little confused by this 'end of the chapter' exercise question from Chapter 2 of "Understanding Digital signal processing" 3rd edition by Richard G. Lyons.

  Consider a continuous time-domain sine wave defined by

  $$x(t)=\cos(4000\pi t)$$

  that was sampled to produce the discrete sine wave sequence defined by

  $$x(n)=\cos(n \pi/2)$$

  What is the sample rate ($f_s$ measured in $\textrm{Hz}$) that would result in the sequence $x(n)$?

There are a few things that I have not understood about this question:

Why is the argument of the cosine function so different after sampling?
I am not sure why the argument of the cosine after sampling does not include $t$ so that we can get $nT$ where $T$ is the sample time.
What determined the argument of the cosine after sampling?
Can I determine the frequency of the sampled signal from the argument of the cosine in the $x(n)$ equation?
And of course a cheeky one, the answer to the book question please :)


texts: What is the difference between multiplying a delta and a step versus convolving a delta and a step?

texts: Seems both will produce another step. there is no difference? Thanks


texts: Does the OFDM MMSE Channel Estimator Require Prior Channel Knowledge?

texts:
The OFDM MMSE channel estimator seems to require a-priori knowledge of the channel to calculate the MMSE estimate.  See equation 8 here: http://cache.freescale.com/files/dsp/doc/app_note/AN3059.pdf.  The autocovariance matrix explicitly includes the actual channel H or g.

However, you presumably don't have this information, otherwise why would you be estimating it?
If you are able to deal with the complexity of the MMSE estimator, how do you implement it in practice without prior knowledge of what the channel is?


texts: How can I lowpass filter by only reducing peak data?

texts:
I have a 2D image, which I want to lowpass filter, with these constraints / quality metrics:

I can not "add" light to the image, so each pixel in the result
should be <= the corresponding pixel in the input.
The lowpass cutoff frequency should be a parameter, to experiment with
Applying this filter repeatedly should not change the result in a significant way.
The time it takes to run this algorithm (5 minutes for a 5MPix image seems reasonable)
Minimizing the amount of light that is filtered away.

Below are some approaches I have tried, together with their shortcomings:

Gaussian filter like normal, then pull the result down to comply with constraint 1. This complies very well with the first 3 points, but reduces a lot more light than necessary.
Fitting "upwards" parabolas through the "low" points and "downwards" parabolas between them to smooth out. This works great in 1D, but applying it first horizontally, then vertically produces bad results in 2D. It takes a lot longer, but not too long for my application. However, repeatably applying this filter will drastically change the result. If the (1D) input is a perfect "downward" parabola (which should not be filtered at all), it will be replaced by 2 "upward" parabolas sitting at the start/end.
Using some other form of 2D "basic" functions and linear solving to find the optimal parameters. This is an idea only currently, not implemented / tested yet.

My domain of experience in signal processing is almost exclusively image processing, so I hope to find alternatives to this problem with the input of experts active in other areas of signal processing.

update 2011/08/18

Based on the current reactions, I decided to make things a bit more clear by adding graphs of a typical input and the results of the 3 approaches I described originally + the suggestions I received so far.
For easy comparison, I used just 1D filtering in these examples.

Input data:

Gaussian filter + bring it down to comply with requirement (1).
You can see that bringing it down results in unnecessary light reductions on the right hand side.

Parabolas
As far as I am concerned, this is pretty much excellent, sadly it does not translate perfectly into 2D by applying first horizontal, then vertical. In this case, you also see I can evaluate the fitted parabolas in floating point resolution, which is a small benefit, but not absolutely required.

Grayscale erosion
Based on the suggestion from rwong, I tried grayscale erosion. I used a structuring element with the same parabolic shape as my "fitted" parabolas. The result is almost exactly the same, so this looks promising. However, there are still a few problems:
 1. My structuring element was not "big enough" (although it was already 801 pixels wide)
 1. I only have "upwards" parabolas, no "downward parabolas to smooth the transition from one parabola to the next.

Median filtering
Only included for completeness, it is not really what I want.

raw data
I pasted the raw input data + the various python commands onto pastebin, so you can experiment with the same data too.
http://pastebin.com/ASnJ9M0p


texts: Slew Rate Filter With Different Smoothing Factors in Each Direction

texts:
I'm smoothing the output of an ADC with a simple first order low pass filter.  Of course the smoothing rate is the same whether the signal is rising or falling.

Preferably the smoothing rate needs to be different when the signal is rising vs when it is falling.  Can anyone refer me to an filter/algorithm that would allow this?

For example if the input signal drops fast the filter output would be smoothed a small amount and drop fast.  When the signal rises (even if it's fast) the output of the filter should rise slowly.


texts: Assumptions for Hurst exponent calculation

texts:
Are there any general assumptions for the calculation of the Hurst exponent?

Does the signal need to be stationary, for example?
Does it depend on the method?
What about the length of the time series, is longer better?

I am interested in using implementing it for the analysis of an EEG signal. I have read many articles but can't find answers for these questions.


texts:
This can be very general question, but I am trying to under the basics. I want to know the difference between the filters which require a physical model in order to perform the filtering  and the ones which doesn't. I read someone one such type of filter is Luneberg. But I could not find anywhere about them. If anyone could explain difference between these two for understanding purpose, it'd be great.

PS- this is not a homework. I am just trying to understand.


texts: Calculating an image's fourier spectrum by hand?

texts:
Suppose I have a $4x4$ image with the following values as its grey-level intensity for each pixel like this:

I want to get its Fourier spectrum. Usually, I would just punch into Matlab and run a fft for everything to be calculated automatically. Now, if I don't have Matlab, and I have to do it by hand and with just an ordinary scientific calculator, can I still do it?

I know I could do the brute force way on the Fourier transform's equation like this:

In this case, $N=M=4$.

So, $f(x,y)$ returns the intensity value of the above image at position $x$ and $y$. Eg, $f(1,1)=1$, $f(2, 0)=3$

But this will be insane to work out the summations so many times on paper. Most calculators can't do summations with complex numbers too.

I have tried doing it with Fourier transform's separability property, but doesn't seem to work as I still have to go through the tedious summations.

What other ways can I use to calculate Fourier transform on a simple image without using Matlab, particularly for the values in the first row of the example image above?


texts: Transfer function between voltage and decibel output for in-ear headphones

texts:
What would be a decent formula for approximating the connection between voltage and decibel output for a typical set of unpowered in-ear headphones?

I'm generating sine waves from the computer: an array of floats between -1.0 and +1.0 gets sent out through the speaker:

X_out = h * sin( f * 2π * t ); // h is the amplitude

Let's say we assume the user to ramp up their volume so that h=1 corresponds to a loud  but comfortable listening level.  So:

h=1 ~ dB_out=80dB

Now, how can I mathematically connect h and dB_out?

I've read that:

v_out ~ h

Somewhere I have seen the formula:

SPL = 20 * Log10( v_out / v_ref )

So, rearranging:

v_out = v_ref * 10 ^ ( SPL / 20 )
k * h = v_ref * 10 ^ ( SPL / 20 )
h = k' * 10 ^ (SPL / 20 )
  = A * exp( SPL ) + B,  for suitable A, B

Is this working valid?

Here is the alternative formulation of the question:

Given a frequency f, and a desired decibel level,  what h do I need to generate that decibel level?


texts: Adding zeroes in between the samples of a discrete frequency domain leads to both zero padding and scaling

texts:
I won’t show ALL the calculations, but I’ll describe the problem fairly enough to understand.
We’ve been taught that upsampling on one domain leads to padding zeroes or added time period which is zero on the second domain.
That means, that if I add zeroes as samples in between the values of a frequency domain (of some periodic  signal) then I’ll get padding zeroes in between the periods.
But, in one class we’ve seen also that adding zero samples in between the values of frequency domain of a continuous time Fourier transform, we can also get the following from calculations:

Which basically means, that if we add zeroes as samples in between the existing samples at the discrete frequency of a periodic signal, then we can both reverse transform it to a scaled version in time domain and/or padded with zeroes version in the time domain. Is that true?


texts: Additive White Gaussian Noise (AWGN) and Undecimated DWT

texts:
One of the benefits of DWT is that it is an orthonormal transform.

There are statements that the energy of noise component mainly concentrates on the high-frequency (detail) part and distributes homogeneously. The energy of noise component is included in more wavelet coefficients with smaller amplitudes, while the energy of the useful signal concentrates on fewer wavelet coefficients with bigger amplitudes. See Wavelet Denoising of Well Logs and its Geological Performance.

As far as I understand, that is due to the orthogonality property of DWT.

Does this behavior of white Gaussian noise also hold true:

for undecimated DWT (like SWT or MODWT).
for less redundant forms of DWT (comparing to undecimated version) - DTCWT, for example.

As far as I understand both of them aren't orthogonal transforms.


texts:
I have written a vocal pitch shifting algorithm which uses LPC to separate the signal into its spectral envelope (predicted by the coefficients) and the pitch information (the residual error).

My understanding (and intuition) was that the residual signal would be unintelligble, but if I listen to the residual created by algorithm it is easily identifiable as a very low volume, slightly distorted version of the original signal.

Is this correct, or should I be getting what I originally expected?

Cheers,
Ben


texts: How to design a digital filter with a frequency dependent phase delay?

texts:
I want to design a filter with a custom phase delay related to frequency. As frequency increases, phase delay should increase.

The time delay as a function of frequency can be expressed as:
$t_d = \frac{L}{C_{ph}} - \frac{1}{f}$

$L$ and $C_{ph}$ are constants of $0.0017$ and $2628$ respectively.

The range of $f$ is $500\mathrm{kHz}$ to $1 \mathrm{MHz}$, the sampling frequency, $f_s$ is $78.39\mathrm{MHz}$.

Thus, when the frequency is $500 \mathrm{kHz}$, the delay should be $1.37\mathrm{\mu s}$ or $107$ samples.

For the pass band, $f$, I have used a filter response of $e^{-i2\pi fd}$, where $d$ is the delay in samples required.

I have designed stop bands with a filter response of zero between $0 \mathrm{Hz}-100\mathrm{kHz}$ and $1.4 \mathrm{MHz}-f_s$.

In MATLAB my code looks like this:

n = 50;
fs = double(fs);
L = double(L);

Cph = 2628;

f1 = linspace(500e3/fs,1e6/fs,100);
f1d = L/Cph - (1./(f1*fs));
%Our array is backwards though innit
f1d = f1d * -1;
f1dz = f1d * fs;
h1 = exp(-1i*2*pi*f1.*f1dz);

fstop1 = linspace(0,100e3/fs, 10);
hstop1 = zeros(size(fstop1));

fstop2 = linspace(1.4e6/fs, 1, 10);
hstop2 = zeros(size(fstop2));

d=fdesign.arbmagnphase('n,b,f,h', n, 3, fstop1, hstop1, f1, h1, fstop2, hstop2);
%d=fdesign.arbmagnphase('n,b,f,h', n, 1, f1, h1);
D = design(d,'equiripple');

fvtool(D,'Analysis','phasedelay');

This is what I get:

Markers are at 500k and 1M.

What am I doing wrong?