texts: Most efficient algorithm for comparing two paths

texts:
I'm working with a Continuous Hidden Markov Model, and using the Viterbi Decoder which finds the best possible path for a given sequence. Here is how my algorithm works:

1) Process the (input) .wav speech signal
2) Extract the MFCC Coefficients (13)
3) Use these values to train the HMM
4) Execute the Viterbi Decoder

And an example input, after step 4 gives the best path of:

0 0 0 0 0 0 0 0 0 0 0 0 0 2 12 6 6 2 2 4 2 2 2 2 2 6 4 2 2 6 6 2 2 2 2 2 2 2 2
8 4 4 2 4 2 2 2 2 6 2 2 4 2 6 2 2 2 2 2 6 2 10 2 6 4 2 2 4 2 2 2 2 2 4 4 2 2 2
6 10 2 2 4 2 2 2 2 2 2 4 2 2 4 2 2 2 2 2 2 4 2 2 2 2 6 2 2 2 2 2 2 4 4 2 2 2 2
2 2 2 2 2 2 2 2 4 2 12 2 6 2 2 2 2 2 2 2 2 2 2 2 2 2 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 2 1 2 2 2 2 2 2 2 6 2 2 6 2 2 2
2 12 2 2 2 13 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 8 4 4 2 2 2 2 2 2 4 2 2 4
2 6 2 6 6 12 2 2 4 2 12 4 2 2 2 10 4 2 6 8 8 2 2 2 2 12 4 2 2 12 6 2 2 2 8 2 4
2 4 2 2 2 12 2 2 2 2 4 4 6 2 2 6 6 2 4 4 2 2 2 4 2 4 10 10 8 2 2 2 2 6 4 2 4 2
2 2 2 4 2 2 4 2 2 8 2 2 2 2 2 2 8 2 6 2 4 2 10 8 4 2 2 8 2 2 4 4 2 2 8 4 2 6 8
4 2 2 12 2 6 2 4 8 4 2 2 2 2 4 2 2 2 2 2 2 2 2 2 2 4 2 4 2 2 2 2 2 2 2 2 2 2 2
8 2 2 2 2 4 2 2 2 2 8 2 2 4 4 2 2 2 6 2 4 2 2 2 12 2 2 2 2 4 2 12 2 2 2 2 2 2
2 8 2 2 2 6 12 2 2 2 2 2 2 2 2

I find the best possible path for each of my training data (MFCC) against the model, carried out in the steps above. This gives an ok result when comparing myself, in that, if I enter a speech file "Yes" as input, then, the training for "Yes" will show a very simular path, whereas, if I entered "Yes" as input and compared against training for "No" then the best path will be different.

My question is this: Which would give me the most efficient way to compare the viterbi paths? I have only researched the Euclidean Distance algorithm, as well as the Hamming Distance so would prefer to use these. I know it's a very subjective question, but, I am asking in terms of being able to match two possible paths.

Hope someone can help me.


texts:
as you can see on page 314 the derivation of the ln[f etc.] equation is done via the matched filter that i have also provided below. i would like to know how this equation is used to derive the ln[f etc.] from its previous probability functions, and how it then finally becomes in the form of ||a-b||<||a-b||

any help would be greatly appreciated. thanks


texts: Determining Parameters of S Function for Image Global Contrast Enhancement

texts:
I know that S-function (membership function from fuzzy set theory) is a commonly used function for image enhancement. I am reading a paper titled - A Novel Fuzzy Image Enhancement Using S-Shaped Membership Function. The S-shaped membership function is defined as follows:

$$\begin{eqnarray*}
S(x;a,b,c)=\begin{cases}
0, & x\leq a\\
\frac{(x-a)^2}{(b-a)(c-a)}, & a<x\leq b\\
1-\frac{(x-a)^2}{(c-b)(c-a)} & b\leq x\leq c\\
1, & x>c.
\end{cases}
\end{eqnarray*}$$

where the parameters are to be found. In the above mentioned paper, it given that the values of a, b, c are entropy(I), median(I) and (max(I)+mean(I))/2 respectively. But there's no explanation for how and why these values were chosen. I tried to look for it in many papers referred by the above paper but they are of no use. Can someone give me some idea please? Thank you in advance.


texts: ISI Nyquist criterion and spectral efficiency

texts:
While reading Bernard Sklar's Digital Communications book I've got confused with ISI Nyquist criterion (not Nyquist rate criterion). There is stated (Sklar, 3.3 section) for symbol rate $R_s$ minimum bandwidth $W = R_s/2$ and for practical cases we have equation:

$W = 1/2 \cdot (R_s + r \cdot R_s)$,

where $r = 0...1$. So theoretically we have 2 symbols per Hz. For BPSK it is 2 bits/s per Hz. For QAM-64 it is 12 bits/s per Hz. So this statement I can't understand, because through my practice I used to think if we use BPSK and we want to achive e.g. 6 Mbps we need to allocate about 6 MHz of bandwidth (actually more because of rounding and synchronization/code redundancy). Where is my mistake in understanding of Nyquist theory?

Upd. Here is the picture from the book, where Nyquist channel is pictured

It seems the bandwidth of such a channel is calculated from $0$ to $1/2T$ and negative frequencies are neglected. Why in the baseband case we need to determine bandwidth in such a way? This can be an answer for the question above.


texts: Variance of Integral of a real white Gaussian Noise Process

texts: In this question, is the answer not equal to infinity ? Answer is mentioned as 6. But my doubt is cant we think of it like a linear combination of many independent random variables each having infinite variance, so the resulting random variable also has infinite variance. This was a question asked in GATE  exam conducted in India for Electronics and Communication stream


texts: Audio texture analysis for visualization

texts:
I'm making a complex music visualizer as a personal project. I've already determined that I will be using a series of modified Goertzel algorithms for frequency detection to place the frequencies on a spiraling piano scale without losing all detail in the bass or calculating unnecessary harmonics.

I'm now looking for analysis for the post-processing effects. I'm looking to estimate the warmth/brilliance (or liveness), and 'grittiness' of the sound as I go. A good example of what I mean by grit is "Dark, Darker, Yet Darker."

How should I approach this problem? I don't have a large, supervised dataset and I only have access to a laptop at the moment, so training a neural net is pretty much out of the question for me. For the record, I'm using Processing (java-based) as my coding environment.


texts: pyrtlsdr: Why is my output different than universal radio hacker?

texts:
I'm trying to build a script to do some signal processing from my RTLSDR. But, the bits I'm getting is not what I would expect.

The preamble using 'Universal Radio Hacker' looks like this:

But after plotting all values my script is seeing (-1 until +1) it looks like this:

Often I'm also getting plots that don't make any sense. At least the last one has some resemblance of the URH preamble, but also a lot of plots don't.

What am I missing.. 'Universal radio hacker' also uses pyrtlsdr, so I know it can work.
- Complex, real, imaginary - Can't say that I understand it fully. Something wrong there ?
- Noise, gain.. Not 100% sure what to do there. URH seems to set gain to 25, but when I do that I don't get anything, so I left it at 'auto'
- I should be getting (almost) all of the samples. In the beginning I did lose some, but after I moved the postprocessing logic into separate processes it seems fine.

My script:

from rtlsdr import RtlSdr
from datetime import datetime
import asyncio
import matplotlib.pyplot as plt
from multiprocessing import Process, Queue

async def get_samples(sdr, q1):
    counter = 0
    timestamp1 = datetime.now()
    # Get this many samples every time
    async for samples in sdr.stream():
        q1.put(samples)
        counter += 1
        timestamp2 = datetime.now()
        if (timestamp2 - timestamp1).total_seconds() > 1:
            # To see if we missed any samples. Should be close to srate samples p/s
            print(datetime.now().strftime("%Y-%m-%d %H:%M:%S"), 'GETSAMPLES -', int((counter * 131072) / (timestamp2 - timestamp1).total_seconds()), 'samples p/s')
            counter = 0
            timestamp1 = datetime.now()

def check_samples(q1, q2, noise, samplesperbit):
    while True:
        values = []
        withinnoisecounter = 0
        signal = False
        if q1.qsize() > 10:
            print(datetime.now().strftime("%Y-%m-%d %H:%M:%S"), 'CHECKSAMPLES - Warning: Queue depth is', q1.qsize())
        samples = q1.get()
        for sample in samples:
            if not signal:
                # Start when signal breaks out of noise range
                if sample.real > noise or sample.real < -noise:
                    print(datetime.now().strftime("%Y-%m-%d %H:%M:%S"), 'CHECKSAMPLES - signal detected')
                    signal = True
                    values.append(sample.real)
            else:
                values.append(sample.real)
                if not (sample.real > noise or sample.real < -noise):
                    withinnoisecounter += 1
                elif withinnoisecounter > 0:
                    withinnoisecounter = 0
                if withinnoisecounter > samplesperbit * 3:
                    # Break when signal is within noise range 3 * samplesperbit
                    break
        if signal:
            q2.put(values)

def decode_fsk(q2, srate):
    # ToDo, now it just plots
    while True:
        values = q2.get()
        plt.figure(figsize=(20, 5))
        # Preamble should be in here
        plt.plot(values[:2000])
        plt.show()

def main():
    noise = 0.9
    # 38400 * 30 to get exactly 30 samples per bit.
    srate = 1152000
    samplesperbit = 1000000 / 38400 / (1000000 / srate)
    sdr = RtlSdr()
    # Just like in URH
    sdr.freq_correction = 1
    sdr.sample_rate = srate
    sdr.center_freq = 868.200e6
    sdr.gain = 'auto'
    # Run check_samples in another thread to make sure we don't miss any samples
    q1 = Queue()
    q2 = Queue()
    p1 = Process(target=check_samples, args=(q1, q2, noise, samplesperbit))
    p1.start()
    # Run decode_fsk in another thread to make sure we don't miss any samples
    p2 = Process(target=decode_fsk, args=(q2, srate))
    p2.start()
    # This is the main loop
    loop = asyncio.get_event_loop()
    loop.run_until_complete(get_samples(sdr, q1))

if __name__ == "__main__":
    main()

Thanks!


texts: Filter - Spatial Domain Versus Frequency Domain

texts:
Note: This is a followup question of

  How can I create a custom band-pass filter?.
  Can you explain the following python source code?.
  Bandpass filter bank.

Suppose, I have an equation of a filter as follows,

$$
h(x, y) = \frac{1}{1 + 0.414 (x + y)^{2n}}
$$

I know that this is an equation in the spatial domain.

So, this is how I obtain a filtered image:

I create a, say, $3\times 3$ kernel using that equation where ($x$: from $0$ to $2$) and ($y$: from $0$ to $2$).
I pad the kernel and the image to the same size.
I apply FFT to kernel and the image.
Multiply them element by element and apply I-FFT to the result.
Convert the result to an image format like Bmp and so on.

.

Now, suppose, I have the following equation,

$$H(u, v) = \frac{1}{1 + 0.414(u + v)^{2n}}$$

I know that this is an equation in the frequency domain.

So, how am I going to create the kernel?

In the previous case, ($x$: $0$ to $2$) and ($y$: $0$ to $2$).

What would be the starting and ending values of $u$ and $v$ in the loop?
What would be the procedure to apply the kernel to an image?
How should I multiply the kernel and the image?


texts:
What techniques may "State of the art iterative numerical solver for optimally smooth gain reduction" refer to?
As given in a description for a particular limiter.
But particularly, I've been wondering, what kinds of numerical approach to limiting it refers to?


texts: What's the shortest possible code that correct all the single bit errors?

texts:
Suppose the original code is a 2-bit binary coding with 3 possible alphabet.
\begin{align}
a&\rightarrow 00\\
b&\rightarrow 10\\
c&\rightarrow 11
\end{align}
Is that $\operatorname{Hamming}(7,4)$ coding? Isn't 7-bit too long?


texts:
I'm doing research on block noise caused by video compression (especially in newer codecs) and I've read the IEEE article about the HEVC deblocking filter (here).

However, when reading the algorithm I'm confused because the formula seems wrong. The first formula to determine if the deblocking filter will affect a block is:

$ |p_2,_0 - 2p_1,_0+p_0,_0|+|p_2,_3 - 2p_1,_3+p_0,_3|+|q_2,_0 - 2q_1,_0+q_0,_0|+|q_2,_3 - 2q_1,_3+q_0,_3|>\beta \ (1)$

It is also written (and shown in a figure) that $\beta $ increases with the quantification parameter QP.

Moreover, later it is written that to choose between normal and strong deblocking, the following equation is used (with $i$ going from $0$ to $3$):

$ |p_2,_i - 2p_1,_i+p_0,_i|+|q_2,_i - 2q_1,_i+q_0,_i|<\beta/8 \ (2) $

Since IEEE is a respectable source and the paper is a couple years old, it would surprise me that there is a typo. However, I do not understand how this formula can work. It does seem to me that the second inequation would always be wrong at least for some $i$ if the first one is true.

Am I misunderstanding how the filter works or is it a typo in the paper?


texts: What are the applications of Non Unique Factorization Domain (Non-UFD) in signal processing?

texts:
Let $m$ be a composite, and $\mathbb{Z}_m$ is the residual ring of modular $m$. The $\mathbb{Z}_m$ is clearly a Non-UFD. What are the applications of the polynomial with coefficients defined on $\mathbb{Z}_m$ in  digital signal processing?

Or, what are the applications of algebraic structures on $\mathbb{Z}_m$ to digital signal processing  and other fields?


texts: how to sample and modulate a gps gold code

texts:
I am trying to simulate the data transmission of a GPS satellite. Telecommunications is not my field so I am struggling a bit with some concepts. I managed to generate the Gold code of the GPS transmitter using shift registers. This creates a pseudorandom 1023-vector with [-1 1] values. Now I want to create a sampled version of this signal up to 5.714Hz assuming the clock of the codes work at 1.023 Mchip/s and then modulate this signal with a sinusoidal carrier of 4.309 MHz sampled at the same frequency. I am not sure how to sample the signal.
From what I've understood, I wrote something like this in Matlab but it's not right :

a = [1 1 1 1 1 1 1 1 1 1 ]
b = [1 1 1 1 1 1 1 1 1 1 ]
for k = 1: 1023
    A(k) = a(10) ;
    B(k) = mod(b(2)+b(6),2) ;
    C_A(k) = mod(A(j)+B(j),2) ;
    Af = mod(a(3)+a(10),2) ;
    Bf = mod(b(2)+b(3)+b(6)+b(8)+b(9)+b(10),2) ;
% shifting
    sg_1 = [ Af a(1:n-1) ] ;
    sg_2 = [ Bf b(1:n-1) ] ;
end
code_vect = (C_A-0.5)*2 ;
% UP TO HERE IT'S OK, I GOT THE C/A CODE, how to modulate it?
code_vect = code_vect(index,:) % this is 1x1023 Gold code vector
sa_fq = 5.714e6     %sampling freq
ca_fq = 4.309e6     %carrier freq
cd_fq = 1.023e6   %code frequency (I am right?)
clk = 1.023e6   %clock rate
time = 1/sa_fq : 1/sa_fq :1023/sa_fq   %time vector
for i = 1:1023
%sampled signal, I'm sure this is wrong but don't know how to do it
code_vect_sa((ceil(sa_fq)*i)-(ceil(sa_fq)-1):(ceil(sa_fq)*i))=code_vect(i)*ones(1,ceil(sa_fq))
end
carr = cos(2*pi*(ca_fq+cd_fq)*time) %carrier
sig_mod = code_vect_sa*carr ; modulated signal

Can you guys give a hint on how to solve it. Thanks in advance


texts: Interpretation of the eigenvalues of the inverse Hessian in a KLT tracker

texts:
I'm a masters student, preparing a seminar in computer vision. Among the topics is the Kanade-Lucas-Tomasi (KLT) tracker, as described in

  J. Shi, C. Tomasi, "Good features to track".  Proceedings CVPR '94.

Here's a web resource that I'm using to understand the KLT tracker. I need some help with the math, as I'm a bit rusty in linear algebra and have no prior experience with computer vision.

In this formula for $\Delta p$ (step 5 in the summary), note the inverse Hessian:

$$\Delta p = H^{-1}\Sigma_x\left[\nabla I \frac{\partial W}{\partial p}\right]^\mathsf{T} \left[T(x) − I(W(x; p))\right]$$

In the article, good features to track are defined as ones where the sum of inverse Hessian matrices have large, similar eigenvalues: $\min(\lambda_1,\lambda_2)>threshold$. I was unable to understand how and where this is derived from, mathematically.

The intuition is that this represents a corner; T get that. What does that have to do with eigenvalues? I expect that if the values of the Hessian are low, there's no change, and it's not a corner. If they're high, it's a corner. Does anyone know how the intuition of cornerness comes into play in the eigenvalues of the inverse Hessian in order to determine $\Delta p$ across iterations of the KLT tracker?

I've been able to find resources claiming that the inverse Hessian correlates to the image covariance matrix. Furthermore, the image covariance indicates the intensity change, and then it makes sense... but I've been unable to find what exactly an image covariance matrix is with respect to an image, and not a vector, or a collection of images.

Also, eigenvalues have meaning in principle component analysis, which is why I get the idea for an image covariance matrix, but I'm not sure how to apply this to the Hessian, as it's usually applied to an image. The Hessian, as far as I understand, is a $2\times 2$ matrix defining the 2nd derivatives for $x$, $y$, and $xy$ at a certain location $(x,y)$.

I would really appreciate help with this, as I've been on it for 3+ days, it's just one small formula and time is running out.


texts: Delay-and-sum beamforming: implementing steering in C

texts:
I am currently trying to understand how delay-and-sum beamforming works.
I have found this website which explains the process of this technique very well in which the following C code is provided:

#include <stdio.h>
#include <math.h>

#define ANGLE_RESOLUTION 500    // Number of angle points to calculate

int main(void)
{
   int numElements = 4;         // Number of array elements
   double spacing = 0.2;        // Element separation in metres
   double freq = 1000.0;        // Signal frequency in Hz
   double speedSound = 343.0;   // m/s
   int a;
   int i;
   // Iterate through arrival angle points
   for (a=0 ; a<ANGLE_RESOLUTION ; a++)
   {
      // Calculate the planewave arrival angle
      double angle = -90 + 180.0 * a / (ANGLE_RESOLUTION-1);
      double angleRad = M_PI * (double) angle / 180;
      double realSum = 0;
      double imagSum = 0;
      // Iterate through array elements
      for (i=0 ; i<numElements ; i++)
      {
         // Calculate element position and wavefront delay
         double position = i * spacing;
         double delay = position * sin(angleRad) / speedSound;

         // Add Wave
         realSum += cos(2.0 * M_PI * freq * delay);
         imagSum += sin(2.0 * M_PI * freq * delay);
      }
      double output = sqrt(realSum * realSum + imagSum * imagSum) / numElements;
      double logOutput = 20 * log10(output);
      if (logOutput < -50) logOutput = -50;
      printf("%d %f %f %f %f\n", a, angle, angleRad, output, logOutput);
   }
   return 0;
}

Beamforming works by steering to the direction that we want to listen to.
Since the above code calculates the delay for -90 to 90 degrees, I thought that it was steering for a sound wave coming from -90 to 90 degrees.

But when going all the way down to that same website, it states "So far the main lobe of the array's sensitivity pattern has been fixed to a single direction; perpendicular to the array."

From that statement, I deduce that the above C code wasn't steering at all? If not, how could I implement it?


texts: Using MATLAB to define a function in signal processing 101

texts:
I want to define a function $y(t)$ that has 2 cases:

$$y(t)=
\begin{cases}
-1 & \text{for} & -T \lt t \leq 0 \\
1 & \text{for} & 0 \lt t\leq T
\end{cases}
$$
where $T$ is an integer.

I already defined $t$ using linspace between $-T$ to $T$ in a resolution of 1000 points.

I was advised to use the command ones, but I dont know how to apply it by the function terms.

Thank you.


texts: Amplitude distribution of Rician faded complex signal

texts:
I have question regarding Rician fading simulation. I am trying to develop Rician fading channel simulator in LabVIEW, but I ran into some problems. I have complex sinusoid as my line-of-signs signal. Then I generated white gaussian noise samples for I and Q components which represent one scatered wave. I added this scattered wave to my LoS sinusoid. According to this I should get a faded signal with folowing amplitude distribution at different signal amplitudes:

When I calculate amplitude and its distribution with histogram I get following result:

I am using the same amplitude values and sigma value. There is significant difference in these results and my questions are:

What could be the problem here?
If I make multiple scatered waves, will it make difference?
Do I need to take into account Rician factor and somehow scale power of generated scatered waves?


texts: Why does the filtered signal amplitude increase exponentially?

texts:
I am new to signal processing. I am having a problem applying a filter to a test signal. I am working with MatLab.
My test signal is a chirp up to Nyquist frequency and my filter is in zeros and poles form. I am using zp2tf  and filter functions but I have amplitudes that increase with frequency up to 3000dB approximately.
The problem is that this amplitudes surely are bad, then I ask for help to understand how can I correct this problem and why it appears.
Next, I paste my code for your review
close all
clear all
clc
fl = 0.05;                            % chirp low frequency
fh = 3.3125;                          % chirp high frequency
fs = 2*fh;                            % signal sampling frequency
dt = 1/fs;                            % delta of time
T = 20;                               % signal durartion
t = [0:dt:T]';                        % time vector
x = 10^(-8)*chirp(t,fl,T,fh);         % test signal
z = [0]';                             % zeros vector (only one)
p = [-3.56047E-01+2.83702E+00*1i,...  % poles vector
     -3.56047E-01-2.83702E+00*1i,...
     -6.28319E-02+0.00000E+00*1i,...
     -8.06237E+00+3.33954E+00*1i,...
     -8.06237E+00-3.33954E+00*1i,...
     -8.06237E+00+3.33954E+00*1i,...
     -8.06237E+00-3.33954E+00*1i,...
     -3.33954E+00+8.06237E+00*1i,...
     -3.33954E+00-8.06237E+00*1i,...
     -3.33954E+00+8.06237E+00*1i,...
     -3.33954E+00-8.06237E+00*1i]';
k = 2.88869E+07*3.57916E+09;          % gain constant
[b,a] = zp2tf(z,p,k);
y = filter(b,a,x);
hfig = figure;
plot(t,db(y));

Thank you for all your help and valuable time!


texts: Filtered Signal vs File Compression Paradox

texts:
1. Original Situation

I have an original signal as a column data matrix n channels data x:mxn (single), with m=120019 the numer of samples and n=15 the number of channels.

Also, i have the filtered signal as a filtered column data matrix x:mxn (single).

The original data is mainly random, centered at zero, from sensor pickups.

Under MATLAB, am using save with no options, butter as highpass filter, and single for casting after filtering.

save essentially apply a GZIP level-3 compression over a binary HDF5 format, hence we could assume the filesize is a good estimator of the information contents, i.e. maximum for a random signal, and close to zero for a constant signal.

Saving the original signal creates a 2MB file,
Saving the filtered signal creates a 5MB file (?!).

2. Question

How it is possible the filtered signal has a bigger size, considering the filtered signal has less information, removed by the filter?

3. Simple Example

A simple example:

n=120019; m=15;t=(0:n-1)';
x=single(randn(n,m));
[b,a]=butter(2,10/200,'high');
 xf=filter(b,a,x);
save('x','x'); save('xf','xf');

creates 6MB files, both for for the original and filtered signal, which is bigger than the previous values due to using pure random data.

In a sense, indicating that the filtered signal is more random than the filtered signal (?!).

4. Evaluative Example

Consider the following:

A filter created from a random signal $x_r$ from gaussian noise $\sim N(0,1)$, and a constant signal $x_c$ equal to $1$.
Disregard the data type, i.e. let's use only double,
Disregard the data sizes, i.e. let's use one column data vector of 1MB, $n=125000$, $m=1$.
Lets consider the $a$ parameter as the Randomness Index for testing: $x=\alpha x_r+(1-\alpha)x_c$, meaning $\alpha=1$ is fully random and $\alpha=0$ fully constant.
Consider a highpass butterworth filter with $w_n=0.5$.

The following code:

%% Data
n=125000;m=1;
t=(0:n-1)';
[hb,ha]=butter(2,0.5,'high');
d=100;
a=logspace(-6,0,d);
xr=randn(n,m);xc=ones(n,m);
b=zeros(d,2);
for i=1:d
    x=a(i)*xr+(1-a(i))*xc;
    xf=filter(hb,ha,x);
    save('x1.mat','x'); save('x2.mat','xf');
    b1=dir('x1.mat'); b2=dir('x2.mat');
    b(i,1)=b1.bytes/1024;
    b(i,2)=b2.bytes/1024;
    i
end
%% Plot
semilogx(a,b);
title('Data Size for Filtered Signals');
legend({'original','filtered'},'location','southeast');
xlabel('Random Index \alpha');
ylabel('FIle Size [kB]');
grid on;

With the following chart as result:

This simulation reproduces the condition of the filtered signal always having a notorious bigger size than the original signal, which contradicts the fact that a filtered signal has less information, removed by the filter.


texts:
I think freqz in a MATLAB toolbox, is the way to obtain DTFT of sequence.
freqz can calculate frequency response of:

H(z)=(Num)/(Den)

We can easily compute the z-transform of any finite sequence x(n) like this:

H(z)=x(0)z^0 + x(1)z^1 + ...

We know in above expression that the Den, is 1.

Recalling that: freeqz(num,den,n) gives the step response in n point.
By x be vector of x(n),

[x1freqz, x1freqzw]=freqz(x,1,3000,'whole');

must gives us  the DTFT.

1)Is it(above statement) correct? what happening if we shift our polynomial?? why?

The second way is to calculate DTFT formula completely, like this:

[X, W]=me_dtft(x1',pi,3000);
figure
title('my')
% plot(W/pi,20*log10(abs(X)));
plot(W/pi,abs(X))
ax = gca;
% ax.YLim = [-40 70];
xlabel('Normalized Frequency (\times\pi rad/sample)')
ylabel('Magnitude (dB)')

function [X, w]=me_dtft(x,whalfrange, nsample)
    w= linspace(-whalfrange,whalfrange,nsample);
    t=0:1:size(x,2)-1;

    X=zeros(1,size(w,2));
    for i=1:1:size(w,2)
        X(i)=x*exp(-t*1i*w(i))';
    end

end

2)I confused, is the range of parameter t in above code, important?

3)Is this implementation correct? Why?

I think there must be dissonance, since the picture:

Telling us something wrong!
These transform are taken from pure sine wave(its code in in the picture), from right you can see the fft , freqz by the manner explained top and(left) the DTFT as explained earlier.

Edit after "Jason R"'s comment:
Ok, also I removed logarithmic scale, since it make me confuse.
After that, intuitively thy are alike, as you can see in next image, but why they are not exactly the same(refer to last image by logarithmic scale?)?

freqz:

[x1freqz, x1freqzw]=freqz(fliplr(XX'),1,3000,'whole');

figure
title('freqz')
% plot((x1freqzw/pi)-1,20*log10(abs(fftshift(x1freqz))))
plot((x1freqzw/pi)-1,abs(fftshift(x1freqz)))
ax = gca;
% ax.YLim = [-40 70];
ax.XTick = -1:.5:2;
xlabel('Normalized Frequency (\times\pi rad/sample)')
ylabel('Magnitude')

Sine sample:

Fs=1000; Ts=1/Fs;
time=0:Ts:1;

Freqs=500;
Xs=zeros(length(Freqs),length(time));

for i=1:length(Freqs)
    Xs(i,:)= cos(2*pi*Freqs(i)*time);
end

XX=Xs;
XX=XX./ max(abs(XX));

figure;plot(time, XX); axis(([0 time(end) -1 1]));
xlabel('Time (sec)'); ylabel('Singal Amp.');
title('A sample audio signal');
sound(XX,Fs)


texts:
I am given the following system and I am checking the additive property:

$$y(t)=x(e^t)$$

where $y(t)$ is the output and $x(t)$ is the input given to the system.

Now this is what I did so far:

\begin{align}
x_1(t)&\rightarrow \text{system} \rightarrow& y_1(t)=&x_1(e^t)\\
x_2(t)&\rightarrow \text{system} \rightarrow& y_2(t)=&x_2(e^t)\\
\end{align}

$$y_1(t) + y_2(t) = x_1(e^t) + x_2(e^t)$$
How will system respond to this input?

$$x_1(t) + x_2(t)\rightarrow system \rightarrow = $$

What I think is that this is what will happen:

$$x_1(t) + x_2(t)\rightarrow system \rightarrow = x_1(e^t)+x_2(e^t)$$

Am I right?


texts:
Not sure if this is the proper place to ask this, but I'm sure someone here will know the answer.

I have this optical camera (TrackIR 4), which, for what you can see on the table, samples images at a rate of 120 frames per second.

Now how would you do a conversion from fps to Hertz to find out the frequency rate at which the camera samples data?

Maybe it's a stupid question, but I wanna make sure I get the right value.


texts: How a diode-detector affects to Frequency Response?

texts:
I have a question related to crystal radio. I can't understand how a crystal detector (diode) shifts the spectrum from the carrier frequency to zero.

The antenna receives bunch of signals and the LC circuit passes the most of the frequencies to the ground excepting the resonance frequency (and the near frequencies). So it's the signal before detector:

The diode passes the positive half-waves to the capacitor/buzzer. So it's the signal after detector:

The capacitor just smooths the high-frequency ripples and the buzzer emits the original demodulated low-frequency signal.
Let's see how it works in the frequency domain.
The antenna receives signals:

The LC circuit extracts the interesting bandwidth:

The non-linear detector as i can understand should shift the bandwidth to zero (and add some other harmonics):

The capacitor/buzzer pair removes the high-frequency harmonics:

The non-linear detector thing i can't understand. How does cutting the half-wave affect to the spectrum? Why the spectrum is shifted down?
I understand how it works in the time-domain. My question is how does it work in the frequency-domain?

I tried to reproduce such behavior in the sound editor. I've generated the sine at 1000 HZ and then modulate it by the 5000 HZ carrier:

So we have 2 peaks at 4kHz (negative frequency shift) and 6kHz. The next step is the detection. So i've reproduce it by using clipping but the result spectrum isn't the same as i expect:

As you can see there is NO peak at 1kHz at all so the low-pass filter is pointless here.
Why here i've got such result? Where am i wrong?
Thanks in advance!


texts: How to find quantization noise in my simulation?

texts:
For college we have to make simple exercise in simulink (matlab). (since my professor is not familiar with simulink, he can't help).

The question itself is, how can I simulate the SNR when uniform quantization is used. To illustrate this:

Because of this quantization there is noise. I need to display this noise. How do i do this?

Additional information: 10 bit quantization. The source has to be 40db.

I know i can get the quantization error when I subtract the output from the input but i do not know how to get the quantization noise.


texts: Gesture recognition from analog sensor value

texts:
Not sure if this is the right place for my question but let me give it a try.
I have a very simple analog sensor (piezo transducer) which is attached to a microcontroller. I'm using this sensor to detect gestures (hand movements).
The mCU sends the sensor data to my PC where I'm running a simple application that plots the data to a graph. In this graph the human eye and brain can recognize easily which gesture is present on the graph, but, of course, the computer has to recognize it as well.

I'm looking for an algorithm or at least a starting point I can detect these gestures. I'm very new to this field so any advice or information would be a great help.

Here you can see a sample screenshot of the graph with two different gestures:

Feel free to set the correct tags for this question.


texts: Convolution in Spatial Domain Is Multiplication in Frequency Domain

texts:
I have to prove convolution in spatial domain is equivalent to multiplication in frequency domain using two matrices.

$$
x(m, n) =
\begin{bmatrix}
1 && 2 \\
3 && 4
\end{bmatrix}
$$

$$
h(m, n) =
\begin{bmatrix}
1 && 0 \\
0 && 0
\end{bmatrix}
$$

When I used matrix method, I got the following result

$$
\begin{bmatrix}
1 && 2 && 0 \\
3 && 4 && 0 \\
0 && 0 && 0
\end{bmatrix}
$$

But when I converted both to DFT using the kernel

$$
\begin{bmatrix}
1 && 1 \\
1 && -1
\end{bmatrix}
$$

and multiply the result, I get

$$
\begin{bmatrix}
1 && 2 \\
3 && 4
\end{bmatrix}
$$

They are not equal. The dimensions are not matching. The first one is $3\times 3$ while the second one is $2\times 2$. Am I doing it correct? Is there something I am missing?


texts: Is it a good idea to whiten a colored noise in order to apply formulas for AWGN?

texts:
I have posted a question before regarding the detection of a signal under colored Gaussian noise:
Energy detection in presence of colored Gaussian noise

Theoretical formulas for setting the detection threshold are widely used for AWGN.

Is it a good idea to just perform a pre-whitening to the colored noise and then apply the well-know formulas for energy detection under white Gaussian noise?

Do you have any happy experience using this approach in practice?


texts:
I am using the Python rtlsdr library to get complex samples from my RTL-SDR USB stick.  I am getting power by squaring the real and imaginary components of these values, and adding them.

What unit is the result in?
If not $\textrm{dBm}$, how can I get to $\textrm{dBm}$?


texts:
I have been going through few research papers around face recognition and I came across two texture extraction algorithms :

Summed-area table by Crow et al. ("Summed-area tables for texture mapping", 1984)
Integral image by Viola et al. (e.g. "Robust Real-time Object Detection", 2002)

Apologizes if I'm wrong, I found that generating texture map in both these algorithms are same.

  i.e., Texture map $I(x,y)$ = Sum of all gray level intensity values
  who lies left and above $i(x,y)$ of the original image.

Please clarify me, what is the main difference between these two algorithms?


texts: Controllability's dependence on B matrix

texts:
Assuming the LTI system:

$\dot{x}=Ax+Bu\\$ , where $x\in R^n$

I understand that

the $B$ matrix is usually given, and it is up to us to derive the optimal $u$ given an initial state of the system and a desired state where we want the system to arrive, in a finite time
a full-diagonal $B$ matrix could guarantee global controllability for the system
any $B$ matrix has to be dependent on the system and the application

In case the $B$ matrix is not a priori given, how could we cook up the "ideal" $B$ matrix, or at least any possible matrix that can render the system controllable?

potential ways, that might or might not be correct:

calculating the controllability matrix for any possible combination and keep those combinations that offer ideally global controllability, or at least the maximum rank of the controllability matrix. But this is computationally too demanding
through Cayley-Hamilton theorem, finding the eigenvalues of $A$ matrix and accordingly deciding for the minimum and maximum number of columns of $B$ matrix. But this implies that we know or can find, the values of the components of the $B$ matrix


texts: How to interpret the effect of different windows in short time fourier transform?

texts: There are many categories of windows, e.g., rectangular, Gaussian, and triangular. What are their effects on STFT?


texts: Rotated object detection by Viola jones Algorithm

texts:
In my project I need to detect rotated and non-rotated faces in real time. I read original articles and try to use common implementation in matlab and openCV. And as I understood Haar features are not rotation invariant. What should I do for successes detection of rotated faces ? I have some ideas about it. Which is more 'real'?

Try to rotate original image with some angle step
Train own classifier on rotated examples
Use any modified Viola Jones Algorithm from research articles
Use another algorithm (like SURF, KLT, Skin based face detection)


texts: Which unit does the power spectrum of microphone output have?

texts: If I sample sound with a microphone and find the absolute square of the (non-normalized) DFT, I get the discrete power spectrum (correct me if I'm wrong). Which unit does it have?


texts:
I'd like to understand when an input signal is time varying or time invariant.
Specifically, I have this function:

$y(n) = x(n) + 0.1 + (0.2)^n + x(n-3)$

I thought that this is time varying because of the $n$ in the exponent. Would anyone mind to clear this up for me? Is this time invariant or time varying?


texts:
var white = Math.random() * 2 - 1;

// Filter to approximate 10 dB per decade
b0 = 0.99886 * b0 + white * 0.0555179;
b1 = 0.99332 * b1 + white * 0.0750759;
b2 = 0.96900 * b2 + white * 0.1538520;
b3 = 0.86650 * b3 + white * 0.3104856;
b4 = 0.55000 * b4 + white * 0.5329522;
b5 = -0.76162 * b5 - white * 0.0168980;
var sample = (b0 + b1 + b2 + b3 + b4 + b5 + b6 + white * 0.5362) * 0.15; // Roughly compensate for gain
b6 = white * 0.115926;

Above is some code that approximates a 10 dB/decade filter, where "white" represents a random number between -1 and 1, producing white noise, while "sample" is the actual value that is passed as a sample to the output buffer.

My question is, how are the coefficients 0.99886, 0.0555179, et cetera actually determined? Including the 0.5362 and 0.115926


texts: Does an asymmetric Bernoulli matrix satisfy the RIP?

texts:
Define an $n\times N$ sensing matrix $A$ by $A_{ij} = 0$ with probability $p$, and $A_{ij} = 1/\sqrt{n}$ with probability $1-p$.  Does $A$ satisfy the restricted isometry property?

For reference, the symmetric case is answered in the following paper:

  R.G. Baraniuk, M.A. Davenport, R.A. DeVore, and M.B. Wakin, "A simple
  proof of the restricted isometry property for random matrices,"
  Constructive Approximation, 28(3) pp. 253-263, December 2008.
  (pdf)


texts: What are the best approaches to characterize an unknown discrete time LTI system? I believe one of the approaches is to input a known input and measure the output to find a transfer function. What are some more approaches to better characterize this black box system?


texts: I want to actuate a motor based on a wav file, but it has too many samples per second

texts:
I am trying to actuate a simple DC motor based on a wav file. At the moment I am using the "cat vs dogs sounds" from Kaggle.

One of the wav files in this DB looks like this:

These audio files are just a few seconds. The sample rate is 16KHz.

My first idea was to just actuate the motor based on the amplitude, but at 16KHz this doesn't make any sense. Also, the signal goes constantly from + to - (this is sort of ok because the motor can move in both directions), but from a human POV, what the motor does doesnt make any sense, because the movements are too rapid. I then tried to "decimate" the signal into a few points per second, but the signals seem to be badly destroyed.

DC motors are actuated using an Adafruit board, using Arduino, and then I have a Python script that interfaces with the Arduino through serial. I have this working fine, and I can make for example a loop from python that increases the speed every few seconds. My idea was to go forward or backwards based on the amplitude of the signal (if it is negative or positive), I am no expert on audio, so this might make no sense. I won't react to other sounds, only to the input wav files.

I don't really know what is the best way of transforming an audio signal of 16KHz into something that can be use to actuate a motor. At the moment I am using a DC motor, but the idea in the near future is to use leds instead, and change their intensity based on the audio signals


texts: Discrete Time Processing of Continuous Time Signals LTI Concept

texts:
I wonder if the input signal (CT) violates Shannon-Nyquist Theorem for a given sampling rate, is there any chance for the overall system not to be LTI although discrete time system is LTI?

Thanks.


texts: Frequency range for set of all discrete-time sinusoids

texts: Consider discrete time sinusoids of the form $$x[n] = \cos(\omega n) \ ,$$ where $n$ is an integer. What frequency range of $\omega$ would constitute all the possible sinusoids? I'm thinking that $0\leq\omega<\pi$ would suffice, but my teacher tells me it's $-\pi\leq\omega<\pi$.


texts:
I am trying to find the impulse response to this system:

At first it looks easy enough, $y(n)=x(n)(-1)^n$, set $x(n)=\delta(n)$, which makes $h(n)=(-1)^n$. But something seems wrong.

What am I missing?

(considering that $(-1)^n=\cos(\pi n)=e^{j\pi n}$, this system seems more complicated than it looks. It seems like $y(n)$ is $x(n)$ with a 180$^\circ$ phase shift.)


texts: How to find PSD (Power spectral density) of spatial data

texts:
I have profiled a surface, measuring the height of peaks and troughs at 0.02mm intervals. I have 1501 data points, the below listed as an example:

$$\begin{align}
x &= 0\,\mu\text m  & y &= 20\,\mu\text m\\
x &= 20\,\mu\text m & y &= 15\,\mu\text m\\
x &= 40\,\mu\text m & y &= 12\,\mu\text m\\
x &= 60\,\mu\text m & y &= 10\,\mu\text m\\
\end{align}$$

How do I generate a spatial as opposed to frequency PSD plot using Matlab? i.e. what function should I use.

Thankyou


texts: Upsampling class activation maps for discriminative feature localisation

texts:
I am currently reading a paper by on learning deep features for discriminative localization where the authors propose to use class activation maps to learn discriminative localised features. The author's replace the layer after the last convolution with a Global Average Pooling (GAP) for the feature maps. They then train the network to learn the weights $w^k_c$ for each GAP feature map $k$ per class $c$.
In the visualisation the activation maps for a particular class, the author's propose to take a weighted sum of the Class activation maps to generate the final class activation feature map. Undoubtedly, the feature maps after the last convolution cannot be used as visualisation because they are smaller than the original input size and the author's propose to upsample those.
My question is what upsampling technique do they use? is it a deconvolution operation to upsample the image? Can simple upsampling really reconstruct the activation maps?


texts: There are many algorithms that can calculate FFT. Even with Cooley–Tukey FFT algorithm, different radix can be used and the algorithms can divided into decimation in time and decimation in frequency. Any comment on how to choose these algorithms in practice?


texts:
This my 1st post here.

I'm using this impulse invariant (MZT) based method to calculate coefficients for biquad digital filter (for RIAA and non-RIAA de-emph.):

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

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

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

Known issue for MZT is magnitude error close Nyqvist at low sample rates - above filter results over -3dB error at 20kHz. I'm trying to find a way to fix this error. When using 96kHz sampling the error at 20kHz is around 0.5dB so the error (range 20Hz-20kHz) gets fixed by rising the sampling frequency high enough.

Not long ago Martin Visanec published an article "Matched Second Order Digital Filters" and he even kindly looked this case and tweaked his method suitable for above filter but the results he posted was not what I'm looking after (though, (without additional pole/zero pair) the result was less than ±0.3dB error in this above case).

Are there fitting methods one could consider to use for to improve the filter at low fs (for filtering system where coefficients are calculated in real time)?

Any thoughts and help with equations would be appreciated.

Juha


texts:
We know the below,

$$
\mathscr{F}\big\{x(t)\big\}=X(f) \tag{1}
$$
$$
\mathscr{F}\big\{x(-t)\big\}=X(-f) \tag{2}
$$
$$
\mathscr{F}\big\{x^*(t)\big\}=X^*(-f) \tag{3}
$$

Now, if for some signal

$$
x(-t)=x^*(t) \tag{4}
$$

Then, is it safe to assume the following?

$$
X(-f)=X^*(-f) \tag{5}
$$

or does it depend on the type of signal?


texts: How exactly does one apply a Hanning Window for a spectral estimate?

texts:
PREAMBLE

I am a graduate student researching evolutionary ecology. My supervisor and I have been trying to learn how to simulate colored noise, discretely, by the $\frac{1}{f^\alpha}$ power law for $\alpha \in \mathbb{R}$, specifically on the interval [0,2].

We found an article that seems to cover a method in great detail.

Link here: Discrete simulation of colored noise and stochastic processes and $1/f^\alpha$ power law noise generation (Kasdin, 1995)

On page 806, we were able to prove to ourselves that by a symmetric autocorrelation function, the formula for the sampled spectrum of an arbitrary segment of a process is indeed,

$$\hat{S}(\omega)\triangleq E\{\tilde{S}(\omega)\}=\int_{-T}^{0}\left(1+\dfrac{\tau}{T}\right) \Big[ \dfrac{1}{T} \int_{t_o+\tau/2}^{t_o+T+\tau/2}  R(t, \tau)\Big] e^{-j\omega\tau}d\tau + \int_{0}^{T}\left(1-\dfrac{\tau}{T}\right) \Big[ \dfrac{1}{T} \int_{t_o+\tau/2}^{t_o+T+\tau/2}R(t, \tau) \Big] e^{-j\omega\tau}d\tau$$

$$=\int_{-T}^{T}\left( 1-\dfrac{|\tau|}{T}\right)R(t,\tau)e^{-j\omega\tau}d\tau$$

Where,

$$R(t,\tau)\triangleq \{x(t+\frac{\tau}{2})x(t-\frac{\tau}{2})$$

Next, the paper computed the spectral estimate for Brownian motion. That estimate in the paper is given as,

$$\hat{S}_B(\omega)=2(1+\dfrac{t_o}{T})\dfrac{1}{\omega^2} - 2(\dfrac{t_o}{T})\dfrac{cos\omega T}{\omega^2} - \dfrac{2}{T\omega^3}sin\omega T$$

Where for $T$ much larger than $t_o$, this reduces to $\frac{2}{\omega^2}$. We were able to reproduce this result.

QUESTION AND PROBLEM

The paper now makes the claim that we have unknowingly been using a rectangular window. We are told we can improve the spectral estimate by using a Hanning window with normalization $8/3$ (After some reading, I've been able to grasp the concept of spectral leakage).

The new spectral estimate for Brownian Motion, with $t_o$ set to 0 is,

$$\hat{S}_B=\dfrac{1}{3}\{  \dfrac{4}{\omega^2} - \dfrac{4sin \omega T}{T \omega^3} + \dfrac{80\pi^4 - 24\pi^2T^2\omega^2 + T^4\omega^4}{T^4(4\pi^2/T^2-\omega^2)^3)} + \dfrac{48\pi^2T\omega sin \omega T - 4T^3\omega^3sin\omega T}{T^4(4\pi^2/T^2-\omega^2)^3)}  \}$$

Try as we might, we have not yet been able to reproduce this result. We tried the following,

$$\int_{-T}^{0}\left(1+\dfrac{\tau}{T}\right) \Big[ \dfrac{1}{T} \int_{t_o+\tau/2}^{t_o+T+\tau/2}  R(t, \tau)\Big] e^{-j\omega\tau}d\tau \cdot HW(x=t)\cdot HW(x=(t+\tau)) + \int_{0}^{T}\left(1-\dfrac{\tau}{T}\right) \Big[ \dfrac{1}{T} \int_{t_o+\tau/2}^{t_o+T+\tau/2}R(t, \tau) \Big] e^{-j\omega\tau}d\tau \cdot HW(x=t)\cdot HW(x=(t+\tau))$$

Where we defined the Hanning Window (HW) as,

$$1-cos((\pi x)/T)^2$$

Done with the following code in Sage:

Q, T, t0, w, tau, t, x = var('Q T t0 w tau t x')
hanning_window = (sqrt(8/3))*(1 - cos((pi*x)/T)^2)
rectangular_window = 1
assume(T>0)

window_1 = hanning_window
window_2 = hanning_window

fneghann = Q*(t+tau/2)*exp(-I*w*tau)*window_1.substitute(x=(t))*window_2.substitute(x=(t+tau))
fposhann = Q*(t-tau/2)*exp(-I*w*tau)*window_1.substitute(x=(t))*window_2.substitute(x=(t+tau))

result = (1/T)*((fneghann.integrate(t,t0-tau/2,t0+T+tau/2)).integrate(tau,-T,0) + (fposhann.integrate(t,t0+tau/2,t0+T-tau/2)).integrate(tau,0,T))
f = (result).expand().full_simplify()
view(f)

Computing the limit as $T\rightarrow \infty$ and setting $t_o = 0$ DOES produce the spectral estimate of $\frac{1}{\omega^2}$ as indicated in the linked paper, but we got something dissimilar to the initial spectral estimate submitted in the article, leading us to believe the result of our limit was luck. I've chosen to not post the result, $f$, because of it's size.

Where did we make a mistake, if any? Did we apply the Hanning Window incorrectly? Is there a mistake in our code? Your help is greatly appreciated.


texts: What is the maximum frequency that can be captured by a given sampling rate?

texts:
Let's say that we want to digitally sample an analog signal and have a sampling frequency of 8000 samples/second. Usually the maximum frequency that can be captured is the half of the sampling frequency; that is the maximum frequency that can be captured is approximately 4KHz for the example above. But what is the mathematical relationship between the two?
Thank you.


texts:
I don't know how to determine line rate and frame rate of the following system:

Height of monitor:            40cm
Screen Aspect Ratio:          16:9
Pixel Aspect Ratio:           4:3
Progressive scan
Viewing distance from screen: 0.6m
Active line period:           15μs
Horizontal Blanking:          5μs
Vertical blanking lines:      20


texts:
I'm trying to create an application using python that is capable of recording an audio signal and detecting short glitches in the signal. For example, I could be recording a 1000Hz sine wave and looking for instances when the audio cuts out for a sample or two. I've been achieving this by removing all the 1000Hz frequencies using a notch filter and searching where the remaining peaks occur. My program seems to work fine when I begin recording a signal, finish recording and then filter the entire recorded signal, search for peaks etc but I would like to get it working in real time, where I can begin recording, analyze 5 seconds worth of audio then move on to the next 5 seconds of audio and so on until I finish recording. Unfortunately my notch filter that worked previously no longer seems to work correctly when I do this.
The code for my notch filter is shown below (I am using scipy for signal processing and pyaudio for recording audio))

f0 = 1000.0  # Frequency to be removed from signal (Hz)
Q = 30.0  # Quality factor
w0 = f0/(fs/2)  # Normalized Frequency
# Design notch filter
b, a = signal.iirnotch(w0, Q)
zi = signal.lfilter_zi(b, a)
z, _ = signal.lfilter(b, a, xn, zi=zi*xn[0])
z2, _ = signal.lfilter(b, a, z, zi=zi*z[0])
y = signal.filtfilt(b, a, xn)

The graph shown below is a filtered 5 seconds of audio. The large peak just after 1 seconds is a genuine glitch but the smaller peaks just after 2,3 and 4 seconds should not be there.

These smaller peaks also continue to occur when I analyze the next 5 seconds of audio, despite there being no glitches at all in this chunk of audio, shown below.

When I begin recording audio, finish recording and then check for glitches, these miniature peaks do not occur at all, they only occur when I try to divde the incoming audio into chunks and check for glitches as the audio is being recorded. Can anyone explain what they are (something to do with harmonics?) and whether or not I can get rid of them? Ideally I would also like to get rid of the peaks that occur at the beginning and the end of every 5 second chunk as well. I dont have much of a background in DSP and I'm still learning a lot about it so I would appreciate any help.


texts: Can anyone tell me how does bit reversal order technique works for FFT when you are trying to break down the FFT to small sizes. Like I want to only for the when the index is odd because that involve some reverse carry propagation, but when the index is even you just add N/2 to it to get your reverse index r. An example of N=8 would be perfect if someone can explain me.


texts: In image processing, how do region growing and clustering differ from each other? There are many papers like clustering based on region growing. I couldn't able to draw a fine line between these two concepts.


texts:
I am struggling trying to understand the 3rd line of the expression below. The second line can be simply visualized by letting the kernel size be 3 and $i=2$ for eg.

Then
$ y_2 = \sum_a x_{2+a-1}k_a = x_2k_1 + x_3k_2 + x_4k_3$

In the 3rd line, there is a change of variable letting $u=i+a-1$. What happens then ? How do I compute $y_2$ now ?

$ y_2 = \sum_u x_uk_{u-i+1} = ?$


texts: What is the difference between "tree-based clustering" and "tree-based classification"?

texts:
Clustering is an unsupervised learning method and classification is a supervised one.

I understand that a tree can be used for image classification since it is based on decision tree in which we have yes/no conditions. we grow a tree to classify a query image whether it belongs to a predefined group or not. But what I cannot understand is clustering using trees. In clustering, how we can adopt a decision tree to group our data in $n$ groups?


texts: Can FFTW perform the FFT on data which is not uniform in time? I can't seem to find a straight-forward answer to this question anywhere online.


texts: Fourier Transforms, Convolution, Cross-correlation: what is their physical unit exactly?

texts:
Let us assume we are talking about real, deterministic, electrical signals $x(t)$ and $y(t)$ (magnitude in Volts).
There are different kind of Fourier Transforms. I made a table to summarize:

NB: By the "<-" symbol, I mean variable substitution.
I tried to figure the physical unit of the output. The $\mathrm{V}\cdot\mathrm{s} = \mathrm{V/Hz}$ for the FT is fine, but I am not satisfied with what I get with the alternative transforms... is that right?
Also, analyzing the units of convolution product of $x(t)$ and $y(t)$:
$$
(x \ast y)(t) = \int\limits_{-\infty}^{+\infty}x(t-u)y(u) \ \mathrm{d}u
$$
... or the cross-correlation of $x(t)$ and $y(t)$:
$$
(x \star y)(d) = \int\limits_{-\infty}^{+\infty}\overline{x(t-d)}y(t) \ \mathrm{d}t
$$
...it would yield $\mathrm{V}^2 \cdot \mathrm{s}$ unit, which is... I don't know... somewhat closer to an energy?
But, as I understand it, the convolution product is actually a signal, typically the output of a filter bank from the original signal.
------------------------------------- EDIT --------------------------------------
In electronics, for the convolution product, the $y$ function should rather be seen as a pattern. Typically  $y$ is an impulse response (often noted $h$), its unit is $s^{-1}$ and then the unit of the convolution product is $V$, which is legit for a signal.
-------------------------------------- END --------------------------------------
On the other hand, the cross-correlation is rather meant to be some kind of inner-products series, so it makes more sense to see it as some kind of energy (the auto-correlation of a signal at $d=0$ multiplied by some coefficient is actually its energy).
------------------------------------- EDIT --------------------------------------
In the cross-correlation product, both $x(t)$ and $y(t)$ are indeed signals.
The above definition is for finite-energy signals, but it changes for finite-power signals:
$$
(x \star y)(d) = \lim_{T\to\infty}\frac{1}{T}\int\limits_{-\frac{T}{2}}^{+\frac{T}{2}}\overline{x(t-d)}y(t) \ \mathrm{d}t
$$
-------------------------------------- END --------------------------------------
So physically speaking, maybe there are some unitary normalization coefficients missing in those formulas?
------------------------------------- EDIT --------------------------------------
For finite-energy deterministic signal $x(t)$, we get:
$$
R_0.E_x = \int\limits_{-\infty}^{+\infty}|x(t)|^2 \ \mathrm{d}t <{+\infty}
$$
For finite-power, we get:
$$
R_0.P_x = \lim_{T\to\infty}\frac{1}{T}\int\limits_{-\frac{T}{2}}^{+\frac{T}{2}}|x(t)|^2 \ \mathrm{d}t <{+\infty}
$$
We may divide those by $R_0 = 1 \Omega$ to get respectively energy or power quantities (hypothetical, relative to 1 Ohm).
Same goes for cross-correlation, depending on whether the signals are finite-energy and/or finite-power.
-------------------------------------- END --------------------------------------


texts: What filter is better with -1/3 2/3 -1/3 or with -1/4 1/2 -1/4 coefficients?

texts: Professor asked me what filter is better -1/3 2/3 -1/3 or -1/4 1/2 -1/4, the answer was second. But what those coefficients means? It was during discussion of linear filter but i'm not sure was it about it.


texts: Distribution of the output of AWGN channel

texts:
$$\mathbf{Y}=\mathbf{X}+\mathbf{W}$$

$\mathbf{X}$ is input, an $n$-dimension equiprobable vector in a finite subset of complex space $\mathcal{A}^n \subset \mathbb{C}^n$ with norm $\sqrt{\rho n}$ (constant power), e.g. BPSK, QPSK, 8-PSK, etc.

$\mathbf{W}$ is circularly complex Gaussian vector $\mathcal{CN}(0, \mathbf{I}_n)$.

I want to know the closed form of $P_{\mathbf{Y}}(\mathbf{y})$. And if possible, the $P_{\mathbf{Y}}(\mathbf{y})$ induced by the QAM modulation of $\mathbf{X}$ (identical average power).

I can numerically calculate
$$P_{\mathbf{Y}}(\mathbf{y}) = \sum_\mathbf{x}P_{\mathbf{Y}|\mathbf{X}}(\mathbf{y}|\mathbf{x})P_{\mathbf{X}}(\mathbf{x}) = \sum_\mathbf{x}P_{\mathbf{Y}|\mathbf{X}}(\mathbf{y}|\mathbf{x})\frac{1}{|\mathcal{A}|^n} = \frac{1}{|\mathcal{A}|^n} \sum_{\mathbf{x} \in \mathcal{A}^n} \frac{1}{\pi^n}e^{-\left\Vert\mathbf{y-x}\right\Vert ^2} \tag{1}$$

where $|\mathcal{A}|$ is the cardinality of set $\mathcal{A}$. But I want to avoid $\sum_{\mathbf{x} \in \mathcal{A}^n}$ because it is cumbersome.

If $\mathcal{A} = \mathbb{C}$ I can use the spherical symmetry trick that
$$P_{\mathbf{Y}}(\mathbf{y}) = \frac{1}{S} \times P_{\left\Vert\mathbf{Y}\right\Vert^2}(\left\Vert\mathbf{y}\right\Vert^2) \tag{2}$$

where $S$ is the surface of n-dimension complex sphere with radius $\sqrt{\rho n}$ and $P_{\left\Vert\mathbf{Y}\right\Vert^2}(\left\Vert\mathbf{y}\right\Vert^2)$ is noncentral $\chi$-squared distribution with $n$ degrees of freedom and noncentrality parameter $n \rho$. Thus I have a closed form of the distribution of $\mathbf{Y}$.

I think this must be a classical problem, it would be great if someone gives me a hint or points out a paper.


texts: Convert low pass continuous time filter design to bandpass, discrete time

texts:
I am trying to convert the continuous time transfer function of a second order lowpass Butterworth filter is given by:

To a bandpass fourth order bandpass digital filter, I first apply the mapping to bandpass:

We therefore obtain the bandpass continuous time transfer function

In order to convert this to a digital filter with with band edge frequencies of 2000Hz and 2819.3 Hz, and a sampling rate of 8 kHz. We first normalise and apply bilinear prewarping

approximating slightly we get

Then applying the bilinear transform

We get the fourth order discrete time filter

However
Checking the results at the normalised angular frequencies 0 and at the bandpass center frequency 3, I find that the frequency response is not as desired. is this an arithmetic mistake or have I missed something, your help would be very much appreciated.


texts: Segmentation of half transparent material, e.g. glass

texts:
I'm totally stuck on an issue regarding the segmentation of glassy objects.
I need to get the object as precise as possible.
My approaches were different.
At first I tried to remove the background, so that only some sharp contours are left.
But that only works for objects which have sharp edges / gradients.
Otherwise the object itself is also removed.
I posted two different images.

I tried to remove the background via morphological operations, like grayscale dilatation and a divison on it.
but it didnt help much.
after it, I tried a k-means with k=3 for getting the modified background separated from the gray and black values of the glass.
That wasn't successfull in some cases, but not overall/in average.
I also tried to make a canny edge detection with an overall blured filter, but that lead to weaker results in form of open contours, a lot of noise, etc. pp.

Canny with automatic threshold results:

testimg = imread('http://i.imgur.com/huQVt.png');
imshow(testimg)
imedges = edge(testimg,'canny');
imshow(imedges);

Same goes for the second image.

As you can see, there is a lot of noise inside and outside and doubled edges from the glas border. Even there are gaps in the edges.

So, I need your advices for getting a general approach for dealing with this problem of half-transparent materials, not for just these two images.

1) Other ideas for removing the background without damaging the object?

2) Other segmentation methods for getting the object separated from the background?

If it's possible, then with Matlab, IPT or statistical toolbox hints. Any other hints are also welcome!

Thank you for your answer in advance.
Sincerely


texts: Should I pass Kalman Filter absolute or offset-from-mean sensor values?

texts:
I'm using Kalman filters to segment the loudness of an acoustic signal from surrounding noise. The problem I've encountered is that muffled or faulty microphones measuring 'silence' (-70dB, -69dB, -71dB, ...) have low variance and so the Kalman filter relies on them as being more stable, "pulling down" the true value.

To solve this, I had the idea to pass in the offset-from-mean of all the signals, which would increase the variance of the "dead sensors" and I could add back the mean to the state estimate delta. But then it occurred to me that maybe the Kalman filter is already doing this by minimising "least squares error"?

A secondary problem I have is dynamically adding or removing sensors without losing my covariance matrix. Currently when sensors leave, I keep the same matrix, but I set it's observation variance R_ij = Infinity, which seems to work. However, I would love to know how to re-arrange my matrix values when adding or removing additional sensors.


texts: Incorporate another dimension to non-negative matrix factorization (NMF)?

texts:
My first data profile has 3 dimensions; atom pattern (I think it is similar to frequency), time, and signal intensity. I can use supervised NMF for source separation straightforwardly. Because the atom pattern can be overlap, the first result is just acceptable.

I can retrieve the second data profile as a Boolean vector. At specific time point, some atom pattern has another underlying pattern which can be describe with a Boolean vector. Are there some good ways to incorporate another dimension of information as a boolean vector into NMF model? I try to find a similar data in audio processing, but not yet found. Any reference or idea are greatly appreciate.


texts: Creating a digital filter, from Laplace to $\mathcal Z$-transform (zero order hold) to code?

texts:
I'm trying to create a digital filter in code(C) but any language is fine. Now I've got an analogue filter that I have represented by an equation in the Laplace domain and I want to try and implement it digitally.

So my filter has this form in the Laplace domain:
$$\frac{as+b}{cs^2+ds}$$

I then use MATLAB's c2d command which uses the zero order hold transformation (I have a really poor grasp on this, so this might be wrong) and it gives me this formula:

$$\frac{\left(5\cdot 10^5\right)z-67}{z^2-z}$$

I tried following an example that I found that used the Tustin's method, though when I use the c2d function in MATLAB with Tustin it gives me an error.

My attempt has been

$$\frac{hz-i}{jz^2-kz}$$

$b_0=-i, b_1=h, b_2=0, a_0=0, a_1=-k, a_2=j$

Then from this I've tried (which is wrong)
\begin{align}
\text{output}&=z_0 b_0+z_1b_1+z_2b_2\\
z_2&=z_1\\
z_1&=z_0\\
z_0&=\text{input}-a_0z_0-a_1z_1-a_2z_2
\end{align}


texts: What is intuitively a syndrome in error correction?

texts: In error correction for example for Reed-Solomon codes there are computed syndromes for received data in decoding process. What is intuitively a syndrome (I don't mean the mathematical expression for it)?


texts:
If I take two analog, band limited signals $ x(t) $ and $ y(t) $ with non-zero $y(t)$ I can define the division:
$$
 z(t) = \frac{x(t)}{y(t)}
$$
In general $z(t)$ may have a very wide bandwidth, also if the two signals had a small bandwidth.
I ideally may lowpass filter out a limited band part of $z(t)$ and make a discrete version $ z[n] $ without aliasing.
If I instead make division from the discrete versions $ x[n] $ and $ y[n] $ (they are band limited so I can also have a sampled discrete version that completely describe them) the result can be very different from the other discrete version $ z[n] $ after the lowpass antialias filter, in fact I have aliasing.
Is there a way to estimate in digital domain the $ z[n] $ that I would obtain making the division in analog domain and filtering out the lowpass part?


texts: Check whether a given equation is linear

texts:
$ a = (x, y) \in \mathbb{Z}^2 $ is given as a pixel. My equation in which $g$ is image(matrix)  is defined as,

$f(x, y) = 56g(x,y)+93g(x−1,y)+92g(x+1, y)−57g(x, y−1)+555g(x, y+1) $

How can we know whether $f(x,y)$ is linear, proving? By saying linear, it is like a filter which is linear. $g(x,y)$ means a pixel at $x$ and $y$. I'm really newbie for image processing and the question will help me a lot what's going on underneath.


texts: The frequency spectrum of a static non-linearity driven by colored noise

texts:
Does there exist any straight-forward theory that can explain what happens to colored noise that passes through a static non-linearity? That is, if you have colored noise
$$v = H u \, , \; u \sim N(0,\sigma)$$
generated by driving an LTI filter $H$ by white noise, and pass it though a static non-linear function $y = f(v)$, e.g. $f(v) = v^2$; what is the spectrum of $y$?


texts:
What is the transfer function of the system described by the following affine input ($x$)-output ($y$) relationship:
$$
y[n] = \alpha x[n] + \beta.
$$
Using the Z-transform we find:
$$
Y[z] = \alpha X[z]+ \beta
$$
and
$$
H[z] = \frac{Y[z]}{X[z]} = \alpha + \frac{\beta}{X[z]}
$$
But I don't know how to take it from here to compute the transfer function! I'm not sure if this can be done. Any help would be appreciated.


texts: I need help in understanding "Nyquist Criterion" definition

texts:
I am researching the split-step parabolic equation and its split step solution as in:
Ozgun, Ozlem & Apaydin, Gokhan & Kuzuoglu, Mustafa & Sevgi, Levent. (2011). PETOOL: MATLAB-based one-way and two-way split-step parabolic equation tool for radiowave propagation over variable terrain.
The idea is to solve equation:

where $F$ is a Fourier transform, $x$ is the horizontal range, $z$ is the vertical range, $k$ is the wavenumber, $n$ the refractive index of air and $p$ is the Fourier transform variable in spectral domain.
The idea is that we have an antenna which propagates in $(x,z)$ directon like shown in the image below

We need to choose the $\Delta z$ that is the altitude increment in a way that avoids alliasing. What I don't understand from the link I pasted is that authors say, that there is a "Nyquist Criterion" such that:
$$z_{max}\times p_{max}=\pi N$$
where $N$ is the Fourier transform size, $z_{max}$ is the maximum vertical height as shown above and $p_{max}$ is the maximum value of the spectral variable. Could anyone try and derive this for me? Tell me where it comes from? Thank you.


texts: Hilbert transform in the frequency domain

texts:
As I know, the Hilbert transform
$$ H(x(t))=\frac{1}{\pi t}\star x(t) $$
in time domain is equal to
$$ -j \operatorname{sgn}(f) \cdot X(f) $$ in frequency domain. so I tried simple example using MATLAB as below,
x=[1,2,7,3];
y1=imag(hilbert(x));
f=[0,1,2,-1];
y2=ifft(-1i*sign(f).*fft(x));

but the result of y1 and y2 are different as below
y1 =
    0.5000   -3.0000   -0.5000    3.0000
y2 =
   0.5000 - 0.7500i  -3.0000 + 0.7500i  -0.5000 - 0.7500i   3.0000 + 0.7500i

just only the real part of y2 is same with y1.
Any one who knows why please explain.


texts:
Can someone tell me if the way I solved this problem is right?

Suppose I have a signal $X_a(t)= 0.5cos(700πt)+0.6cos(720πt)+0.1cos(780πt)$ where $F_s=8000Hz$

I want to determine the minimum samples which must be available in order the frequency resolution of DFT be ok.

I did this:

Firstly calculate each frequency of the signal:
$F_1=350Hz$, $F_2=360Hz$, $F_3=390Hz$
Calculate the differences: $F_2$-$F_1$=$10Hz$ ,$F_3$-$F_2$= $30Hz$ , $F_3$-$F_1$= $40Hz$
$ΔF_{min}$=$10Hz$
$ΔF$ = $ \frac{f_s}{N} $, so $N\geq\frac{8000Hz}{10}$=800 samples

Can someone tell me if my solution is right and explain why we choose the minimum difference

Thanks


texts:
It will be great if someone can explain me what exactly means "real impulse response".
Further more , what is the effect of such a response on ROC (Laplace plane) and in particular if its restrict the possible area of the poles and zeros?
Thanks


texts: How to mix two audio files recorded with different settings to obtain a realistic result

texts:
I am a software engineer, currently working on a signal-processing related project, which I am new to. For evaluating the algorithm(which is written in R), I need to mix (merge) two different audio signals synthetically. The question is how to mix these two signals and obtain a “realistic” result?

To be more specific, consider the following scenario:

I have recorded Signal_A, and downloaded signal_B from the internet (therefore two different microphones/settings are used for recording).
The sound pressure level of both signals are known( I measured the sound pressure level of Signal_A in 1m distance, and also know the sound pressure level of Signal_B in the same distance).
As far as I know, the value recorded/represented in each waveform of any signal is relative to that corresponding signal, and since two signals are recorded using different settings we can expect that one might be more amplified than the other. Therefore a simple mixing of the two won’t give us a realistic result.

Now the question is: is there a way to manipulate recorded signal_A and signal_B programatically, assuming the fact that the SPL values of the original signals at 1m distance is known, and get the merged signal as they were playing simultaneously and recorded at 1m distance? i.e. to calculate factor ‘a’ and ‘b’ in a way that the realistic merged signal is formulated by:

Mereged_Signal = a*Signal_A  + b*Signal_B

Thanks in advance!


texts:
I am doing simulation of FMCW radar. I am applying signal processing chain on raw radar data. Here in the code below, I am taking fft of data and plotting it. Moreover, I am doing zero padding by multiplying the signal length by 1000. I should be getting beat frequency at 2.020 MHz instead I am getting at 2.019MHz. There is something little probelem but I do not know what.

Is it the correct way for Zero padding? Moreover, am I defining frequency axes correctly? Thanks everyone in advance

% second-order butterworth
[b, a] = butter(2, [w1 w2], 'bandpass');

load('fb2030');                                   % loading the data
x = fb2030(30,:);
% filtering
y_filt = filter(b,a,x);                           % filtering the received signal
nfft = length(y_filt)*1000;                       %%%% Zero Padding
res = fft(y_filt,nfft)/ nfft;                     % normalizing the fft

f = fs/2*linspace(0,1,nfft/2+1);                  % choosing correct frequency axes
res = res(1:nfft/2+1);                            % amplitude of fft(taking the half length of nfft)
figure,plot(f,2*abs(res));
xlabel('Frequency in MHz');ylabel('amp')


texts: Why is the Periodogram Used as an Estimate of PSD?

texts:
Why is the periodogram method as given by Schuster (1898) - which is not a consistent estimate of the PSD - still used, as opposed to consistent non-parametric methods for estimating the PSD like Welch?
I'm looking at the Signal Analyzer app in MATLAB which seems to switch between the periodogram method with a Kaiser window and the Welch method (also with a Kaiser window).
From Proakis & Manolakis (2006), pg 970: However, in general the variance of the estimate $P_{xx}(f)$ does not decay to zero as $N\rightarrow\infty$." ... "Hence we conclude that the periodogram is not a consistent estimate of the true power density spectrum (i.e., it does not converge to the true power density spectrum)."


texts: Autocorrelation of the product of deterministic and random signal

texts:
I was wondering how to calculate the autocorrelation of a deterministic signal $x(t)$ multiplied by a stochastic process $M(t)$, whose autocorellation $R_M(\tau)$ is known a priori. In my case, $x(t)$ is a truncated monolateral exponentially decaying function.

I suppose that the result of such multiplication  $y(t) = x(t) \cdot M(t)$ is again a stochastic process, but when approaching the calculation of the autocorrelation of $y(t)$ I obtain something that is not even, therefore I suppose I am making some mistakes. I know that the definition of autocorrelation for deterministic signals is different from the one of stochastic processes, but I do not know how to connect the two of them.


texts: Is acceleration noise modelled differently in EKF and UKF Kalman Filters?

texts:
In a lecture on basic Kalman filter, I came across the following assumption about acceleration noise.

Every component of the noise vector $\nu$ is itself a product of time and acceleration values. Nothing is known about the accelerations $a_x$ and $a_y$ but each component of $\nu$ is assumed to have zero mean and every pair of components $\nu_j$, $\nu_k$ has some known variance $\sigma^2_{\nu_j\nu_k}$. As a consequence the last term in each of the four equations is set to zero (not shown here). The covariance matrix Q contains time though as in the screenshot below.

In a latter lecture on Unscented Kalman Filter the assumption was changed as below.

Here the components of noise vector $\nu$ do not contain any effect of time. These time-free components are assumed to have mean zero and covariance $Q$. As a result their effect on the mean state vector $x$ cannot be set to zero. Also $Q$ is independent of time.

Is it standard practice to make different kind of assumptions?


texts: Why do I keep seeing this equation for $E_b/N_0$

texts:
I keep seeing $$\frac{E_b}{N_0} = \frac{E_\text{s}}{\rho N_0}; ~~ \rho=\log_2{M}$$
But my own calculation is:
$$\frac{E_b}{N_0} = \frac{E_s}{N_0}\frac{1}{k}$$
Where $k = \frac{\rho n}{ms}$, where $\rho = msr$ is the spectral efficiency, $m$ is the modulation efficiency, $s$ is the percentage of carriers that carry data from the actual coded bitstream, $r$ is the coding rate $\frac{k}{n}$, $n = sNm$ is the total number of coded data bits in the symbol and $k$ is the total number of information bits in the symbol and $N$ is the number of subcarriers.
Wikipedia shows the following:

There are multiple issues with this

It uses $\rho$, which is spectral efficiency, and actually refers to it as spectral efficiency but then substitutes it with $\log_{2}M$ which actually is the modulation efficiency
It says that 'this is the energy per bit, not the energy per information bit', to follow up its use of modulation efficiency now instead, but there's a contradictory equation that is using the net bitrate i.e. the information rate $f_b$ in the section above
The contradictory equation, which is correct, produces my formulation, and the formulation that I've seen on a few sources i.e. this one: https://uk.mathworks.com/help/comm/ug/awgn-channel.html

$$\frac{E_b}{N_0} = \frac{\frac{P_C}{f_b}}{\frac{P_N}{B}}  = \frac{P_C}{P_N}\frac{B}{f_b}$$
and
$$\frac{E_s}{N_0} = \frac{\frac{P_C}{F}}{\frac{P_N}{B}} = \frac{P_C}{P_N}\frac{B}{F}$$
Which shows how $E_b/N_0$ differs in that it has a $f_b$ instead of an $F$ term. Therefore $F$ needs to be timesed by something that produces $f_b$, and that thing is $k$. $Fk$, the baud rate times the number of information bits in the symbol, is the information rate.

The equation I keep seeing across multiple sources as well as wikipedia is dividing the energy of a symbol by modulation efficiency rather than the number of bits in the symbol, which does not make semantic sense if the symbol has multiple subcarriers, and if it doesn't, this would only give the energy per gross bitrate transmission bit.

My guess as to what they have done wrong is calling $\rho$ $\log_{2}M$ or using $\log_{2}M$ in these equations at all, because AFAIK $E_b/N_0$ is with respect to the information rate. If they use actual spectral efficiency, then it is correct IF the symbol has one carrier. If it has multiple then they'd have to times the single carrier spectral efficiency $\rho = mr$ by $sN$ to get $\rho sN$ i.e. $k$


texts:
I have the following code:

sigma_th = 0.08;
ph_noise = sigma_th*randn(1,2^22);
[b,a] = butter(1,4e6/1e9,'high');
ph_noise_out = filter(b, a, ph_noise);
figure(1);
[psd1,w] = pwelch(ph_noise,blackmanharris(2^22/8),[],[],100e9,'onesided');
semilogx(w,10*log10(psd1/2));grid on;
hold on;
[psd2,w] = pwelch(ph_noise_out,blackmanharris(2^22/8),[],[],100e9,'onesided');
semilogx(w,10*log10(psd2/2));grid on;
[psd3,w] = pwelch(ph_noise_out - ph_noise,blackmanharris(2^22/8),[],[],100e9,'onesided');
semilogx(w,10*log10(psd3/2));grid on;

Here, I generate a noise signal ph_noise and high pass filter the noise using a butterworth filter to get the output ph_noise_out.
If I check the power spectral density of the ph_noise and ph_noise_out, I see that the output spectrum is filtered at low frequencies but is exactly the same at high frequencies. But when I plot the difference between the two noises, I see that the difference falls down as 20dB/dec. Why is the difference falling so slowly? How can I make it fall faster?
PS. The blue curve is the ph_noise, orange is the filtered noise ph_noise_out and yellow is the power of the difference signal.


texts:
I'm trying to get the spectrum of my signal in dBFS units. I'm not sure what is the refernce that I should devide by. I computed the power spectrum as follows:

1) FFT normalized to the number of samples

p = abs(fft(signal))/(No_of_samps/2);

2) and then I calculated the power spectrum for half of the samples (one-sided)

p = p(1:Nsamps/2).^2;

3) I  tried to represent it in dBFS as follows:

p_dBFs = 10*log(abs(p)/(vref/2)^2);

Where vref is my reference p-p voltage in the ADC.

I'm just not sure whether I should devide by the reference voltage to get values in dBFS, like when I apply a FS signal I would like to see the highest peak in the spectrum at exactly 0 dBFs.

Note: The time-domain signal is also multiplied by a hamming window before computing the FFT.

Any help is appreciated.

Sorry for the long post. Thank you guys very much.


texts:
I have a doubt, it may be silly. Please explain me.

In K-L Transform, For forward transform we have the below equation

For inverse transform we have the below equation

In reverse transform why the 'A' transpose is not common to both y and m(mean of x) as in Forward transform equation?


texts: question related to something in karlin and taylor stochastic processes one text

texts:
This question is essentially a question about something in Karlin and Taylor's Stochastic Processes One text in the spectral chapter. Since this is a DSP list, Karlin and Taylor may not be so popular so I will describe the setting which I think I've seen in other texts also.

They define the random variable $X_n =  A  \cos(\omega n) + B \sin(\omega  n)$.

So, $A$ and $B$ are zero mean normal random variables. $\omega$ is fixed at some frequency between $0$ and $\pi$. The process is on  the integers so the process start at some value of $-k$ and continues through to positive $k$. So $k$ can be any integer and $X_n$ represents a discrete stochastic process.

They then say that the $X_n$ process is stationary in the mean. Evidently this is so obvious that there's no need to explain it. I don't get it. Well, I sort of half get it.

A and B are  mean zero so, okay the mean will always be zero no matter what the values of the cos or sin term.  But since $n$ is increasing by 1 as the process proceeds, isn't the trajectory of $X_n$ changing at each $n$ because it multiplies omega ? So, because A and B are zero, stationary in the mean works. But, say A and B were both 1.0 say. Would the $X_n$ process still be stationary in the mean ? Since $X_n$ and $X_{n+1}$ will be at a different frequencies, $\omega  n$ and $\omega (n+1)$ respectively, doesn't that make the  mean different at times $n$ and $n + 1$ ? so it can't be a stationary process.

So, I guess my question is: For this process,  does stationary in mean only hold because $A$ and $B$ are zero in mean or would it work no matter what their means were ? If the latter is true, then I'm totally not understanding the process. Thanks.


texts: Is the Laplace transform a special case of Fourier transform? (Not the other way around)

texts:
Always had a thought about why Laplace transform reveals the transient properties of the system?
My doubt is based on the following fact,
Fourier transform is given as

\begin{equation}
\mathscr{F}\left\lbrace f(t)\right\rbrace = \int_{-\infty}^\infty f(t) e^{ -j \omega t}  dt
\end{equation}

Where Mathematically and intuitively we believe that the angular frequency $\omega$
takes only real value.

What if, instead of taking real angular frequencies, if the variable
$\omega$
assumes a complex angular frequency in the form
$\beta - j \alpha$
, then,

$$
j \omega t = j (\beta - j  \alpha) t = (\alpha + j \beta ) t = s t
$$

While taking Fourier transform w.r.t $\omega$, the quantity  $\beta$ will be real angular frequency in radians per second and $\alpha$ will be the $\textbf{imaginary angular }$ frequency in radians per second.

\begin{equation}
\int_{-\infty}^\infty f(t) e^{ -j \omega t}  dt = \int_{-\infty}^\infty f(t) e^{ - s t}  dt = \mathscr{L}\left\lbrace f(t)\right\rbrace
\end{equation}

Hence is it mathematically correct to consider bilateral Laplace transform as a special case of Fourier transform (not the other way around) when $\omega$ takes a complex angular form $\beta - j  \alpha$ ?  I believe the fact that $\omega$ can take complex values is the reason why we get transient properties of the system when using Laplace transform.


texts: Better understanding of downsampling (decimation) and upsampling (interpolation)

texts:
Although some questions were asked about this topic, I have not seen any that answers all the basic questions, that is why I took the liberty to ask more about this.

I suggest to limit, (or at lest to start by), 1D signals. Of course if someone wants to extend they are welcome.

For example, in some lectures it is mentioned that

an ideal 1D interpolator is an upsampler (by the wanted factor, e.g. 2) in cascade with an ideal lowpass filter i.e. an sinc.
Decimation would be to opposite. A convolution with a filter (sinc?) then a downsampler (by the wanted factor e.g. 2).

In Matlab, given a sine wave at a given frequency for example, how would you proceed to downsample it to a lower frequency? And could you make the link with aliasing ?

I have started to do the following, "manually" coded:

    %Working with audio frequencies, for example,
    %to listen the results with audioread() function afterwards
Fs=8000;%sampling frequency
t=0:(1/Fs):1;
f=4000;
x=sin(2*pi*f*t);
x_dwnsamp=zeros(length(x));
%decimate x to see what happends wrt aliasing
decimation_step=2;
for i=1:decimation_step:length(x)
   x_dwnsamp(i)=x(i);
end
figure;
subplot(1,2,1);
plot(abs(fft(x)));
subplot(1,2,2);
plot(abs(fft(x_dwnsamp)));

This was an attempt to resample the signal x which is just at the Nyquist freq. (4KHz) fora sampling freq. of 8KHz.


texts: Why is an audio compressed in FLAC(Free Lossless Audio Codec) changing its bit rate?

texts: I just recognized that the bit rate of an audio file compressed in FLAC is changing during play. On the other hand, MP3 file has a fixed bit rate such as 128, 192, 320kbsp etc. But why FLAC is not? Thank you for your time.


texts: Algorithm to change an already sampled guitarstring to sound like it was a bit damped by the hand

texts:
I'm trying to simulate the damping of the hand on a guitarstring on an already recorded/sampled open guitarstring sound. I've been trying to use low pass filter and had a moving frequency range but that didn't make it sound like a palm muted string, just the loss of higher frequencies.

Could someone help me find good material on this, that a human could atleast grasp a bit? It's going to be implemented in C++ and I have been searching and found almoust everything about the karplus-strong string algorithm, but that's not what I want.. I do want the damping part implemented on a sample of an already recorded real played string.

My question was first asked here.. I hope that will not be considered as crosslinking now?? https://stackoverflow.com/questions/19667049/algorithm-to-change-an-already-sampled-guitarstring-to-sound-like-it-was-a-bit-d


texts: Why correlation is not associative in the context of image processing?

texts:
I can see most of the places it says "convolution is associative, while correlation, in general, is not".

Denote $*$ convolution operator, let's say you have an image $f$, which you need to convolve with $g$ and then with $h$ :
$$f∗g∗h=f∗(g∗h)$$

But when I see the difference between convolution and correlation, correlation is equal to convolution but after flipping the kernel or window.

Can somebody explain with an example why correlation is not associative?


texts:
This is more like an optimization problem but any solution is appreciated.

I have a data set with input specifying power(demand) to be generated for a particular time period(TP).

Input:
Time  ---    Power
01 to 02             1000
02 to 03  2000
03 to 04  2000
04 to 05  1500
etc.,

But the power can't just be increased from 1000 to 2000 at 01:59:59.
Here there will be some ramping rate(gradient) which I want to minimise by keeping the imbalance to 0.  Imbalance may arise if you deviate from what is needed, here power cannot be more or less than what is required in that Time Period.

In summary I have to generate the black line(which is output) below (red line is the input).

I am trying to get an output which can be a curve(black line) equation or
the power to generate at the deflection points(say at around 1:45 we generate 800 then at 2:15 it can be 2200 which will keep the imbalance to 0).

When I tried to treat this as an optimization problem it may look something like this:
Objective :   Minimise gradient f(x) = Yj-Yi/Xj-Xi between 2 points (Xi,Yi) and (Xj,Yj)
Subject to constraints:

average of Yi = Power value(redline in that Time Period)

But I am struggling to complete it how to represent mathematically to be able to implement using programming like R's lpsolve/optimize or Python's scipy optimizers.

Any help on the approach to take would be greatly appreciated.


texts: Translating Matlab scrambler function to C

texts:
I am trying to implement the following matlab function in c:
Mathworks scrambler Function Documentation

This is the code I've written:

#define N 5
#define MAX 70
#define BITS 7
//char* msg = "This message is Coming to You!";
//char* msg = "Hello world 000\nHello world 001\nHello world 002\nHello world 003\n";
//char* msg = "Hello world 000\nHello world 001\n";
char* msg = "Hel";

uint8_t current[N-1] = {0};
uint8_t to_Xor[N-1] = {0};

uint8_t bin[MAX*BITS];

int main()
{
    int len=strlen(msg);
    int kk = 0;
    int m = 0;
    uint8_t letter = 0;
    uint8_t chars = 0;

    for(int nn = 0; nn < len*BITS; nn++)
    {
        bin[nn] = 0;
    }

    for(int ii = 0; ii < len; ii++)
    {
        letter = msg[ii];
        for(int jj = 6; jj>=0; jj--)
        {
            bin[kk]=(letter&(1<<jj))>>jj;
            //printf("%d\n",bin[k]);
            kk++;
        }

    }

    uint8_t xor = 0;
    uint8_t polynom = 0x0;
    uint8_t scrmbld_output[len*BITS];

    for(int c = 0; c < len*BITS; c++)
    {
        scrmbld_output[c] = 0;
    }

    for(int bin_ndx = 0; bin_ndx < len*BITS; bin_ndx++)
    {

        for(int n = 0; n < N-1; n++)
        {
            if((1<<n)&polynom)
            {
                xor ^= current[n];
            }
        }

        for (int k = N-2; k > 0; k--)
        {
            current[k]=current[k-1];
        }

        current[0] = xor ^ bin[bin_ndx] ^ 1;

        scrmbld_output[bin_ndx] = current[0];
        xor = 0;

        printf("original: %d scrmbld: %d\n",bin[bin_ndx], scrmbld_output[bin_ndx]);
    }

    int k = 0;
    m = 6;

    for(int i = 0; i < len; i++)
    {
        for(int j = 6; j>=0; j--)
        {
            chars |= bin[k]<<m;
            k++;
            m--;
        }
        printf("%c", chars);
        chars = 0;
        m = 6;
    }

    printf("\nlen: %d", len);
    return 0;
}

The results I get from matlab are different.

If I create the following vector and scrambler in matlab:

a = [ 1     0     0     1     0     0     0     1     1     0 ];
scrmblr = comm.Scrambler(2,[1 1 1 0 1], [0 0 0 0]);
ans = scrmblr(a');
ans' = [1     1     0     0     1     0     1     0     1     1]

but what I'm getting is:

original: 1 scrmbld: 0
original: 0 scrmbld: 1
original: 0 scrmbld: 0
original: 1 scrmbld: 0
original: 0 scrmbld: 0
original: 0 scrmbld: 0
original: 0 scrmbld: 1
original: 1 scrmbld: 1
original: 1 scrmbld: 1
original: 0 scrmbld: 1
original: 0 scrmbld: 0
original: 1 scrmbld: 0
original: 0 scrmbld: 1
original: 1 scrmbld: 0
original: 1 scrmbld: 0
original: 1 scrmbld: 1
original: 0 scrmbld: 1
original: 1 scrmbld: 1
original: 1 scrmbld: 0
original: 0 scrmbld: 1
original: 0 scrmbld: 0
original: 0 scrmbld: 0
original: 0 scrmbld: 0
original: 0 scrmbld: 0

I don't see where I went wrong. Any help greatly appreciated, thank you!


texts: Comparing arithmetic complexity of FFT radix-2 and convolution

texts:
Let's assume we have a discrete linear time invariant system and we have a real signal $x[n]$ with length N=50 as input for the system.

The impulse response $h[n]$ of the system is considered to be also a real signal with length M=10.

We want to calculate the output signal $y[n]$.

There are two ways we can do this: We can either convolute the signals $x[n]$,$h[n]$ in the time domain or use FFT radix-2 and work on the frequency domain.

I want to compare the arithmetic complexity of these two methods in number of real multiplications.I know that the arithmetic complexity of convolution is $N \cdot M$ since $N \cdot M$ real multiplications are required in order to calculate $y[n]$. However, I am bit confused with FFT-radix 2.

I know that it requires $μ(Ν) = \frac{N}{2}\log _{2}N$ complex multiplications. Every complex multiplication is equivelant to 4 real multiplications. So how many real multiplications are required to calculate $y[n]$? Do I have to calculate first both $X[k]$ and $H[k]$ by using $\frac{50}{2}\log _{2}50 + \frac{10}{2}\log _{2}10 $ real multiplications? And if so, how do I go on?


texts: Continuously changing pitch of stable-pitch audio sample

texts:
I have an 44100-rate 4 second sample of a legato violin playing at 440 Hz.  I would like to modulate this sample so that it starts at 440 Hz, then transforms its frequency continuously and exponentially to 660 Hz, while maintaining the original timbre.
I have tried pitch shifting directly from 440 to 660 and this works fine, but trying to simply cut the sample into extremely short bits and increase the pitch shift for each consecutive bit resulted in a lot of artifacts.
Do you have any suggestions?
Bonus points if it would be easy to implement in python!


texts:
this is a basic question, but i barely find any tutorial to explain this.  Most of tutorial dont talk about why use IFFT and FFT, they just say "oh, at transmitter, we just take this information time domain as frequency domain, then we need IFFT". I have no idea why we want to treat as frequency domain when it was time domain. There must be some good reasons, right?

So, why we need the FFT pair for Tx & Rx for OFDM? Why not something else ? Someone could help me with some pictures?  Thank you again, friend.


texts: FIR amplitude response using Hr_type3 function in MATLAB

texts:
I am trying to design a type 3 FIR filter and generate its amplitude response from the function Hr_type3 (mentioned in Digital Signal Processing Using MATLAB by J. Proakis), but getting an
Error  : Unrecognized function or variable 'Hr_type3'.
Can any one help me on this?
I am trying to build a digital differentiator from this code.
M = 21;
alpha = (M-1)/2;
n = 0:M-1;
hd = (cos(pi*(n-alpha)))./(n-alpha);
hd(alpha+1) = 0;
w_ham = (hamming(M))';
h = hd .* w_ham;
[Hr,w,P,L] = Hr_type3(h);


texts:
I have a question regarding discrete inverse Fourier transform, and no answer I found on the internet seem to be satisfying. This might be because I do not fully get some of them, so please excuse my ignorance.

Suppose you have full knowledge of a frequency-domain function $\hat{f}(\omega)$. It is a continuous, well-behaved bump function. You want to compute its inverse Fourier transform $f(t) = \Re \int_{-\infty}^{+\infty} \hat{f}(\omega) e^{j\omega t} d\omega$, but this cannot be done analytically, so you compute a discrete inverse Fourier transform. You know the interval $ \left[ \omega_{\min}, \omega_{\max} \right]$ on which it is reasonable to numerically compute the IDFT, so you now have to sample $\hat{f}(\omega)$ on this interval.

My question is the following: is there an equivalent of the sampling theorem for the inverse Fourier transform ? Is there a generic rule to sample my frequency-domain signal, as there is for time-domain ? I noticed that I obtain something similar to aliasing when I undersample $\hat{f}(\omega)$.

Side note : all of this is used in a finite element code to compute a complex source. As such the IDFT will be performed at least hundreds of thousand of times, so the question to limit its cost.

Thanks in advance !


texts: gnuradio extremely sluggish with QT (but seems OK with WX)

texts:
I am evaluating gnuradio with Ettus USRP X310 and 2xUBX-160 daugterboards, connected via GBit ethernet to a new Lenovo laptop running Win10 and gnuradio 3.7.13.5.

I am testing it with the most basic setup: connecting "RX/TX" back-to-back to "RX2" via straight SMA cable (Remark: I know this is risky and I am careful not to exceed the allowed input power). This simple gnuradio diagram sends out a 10kHz tone at baseband; sample rate is set to 1MS/s and RF carrier 1GHz:

This somewhat works but whenever I change something in the TX path (e.g. the baseband frequency or BB attenuation) it takes way too long for the plots to update: I measured more than 8s(!) from the time I press enter to the moment the plots update.
This does not seem to be the case for the RX side (e.g., when I change the RX gain, the plots update relatively quickly).

I also just tested it with QT and WX: It only occurs with QT; with WX the plots update relatively quickly.

Since WX is "deprecated" I would prefer to use QT.

What could be the cause of this?
Is there any basic example (that does TX and RX at the same time) I could start with to have a "correct" baseline?


texts:
To find the transfer function of a channel we say that it is
$$
H(s) = \frac{y(s)}{x(s)}|x(s)=0 for <0
$$
Why we do not define it like
$$
h(t) = \frac{y(t)}{x(t)}
$$


texts:
I need to use Parseval's Theorem to find
\begin{equation}
S\:=\:\sum _{n=-\infty }^{\infty }\:\left[\left(\frac{\sin\left(\frac{\pi }{4}n\right)}{2\pi n}\right)\left(\frac{\sin\left(\frac{\pi \:}{6}n\right)}{5\pi \:n}\right)\right]
\end{equation}
I used the fact that
\begin{equation}
x_1\left[n\right]\:=\:\frac{\sin\left(\frac{\pi }{4}n\right)}{2\pi n}\:=\:\frac{1}{8}\operatorname{sinc}\left(\frac{n}{4}\right)
\end{equation}
and
\begin{equation}
x_2\left[n\right]\:=\:\frac{\sin\left(\frac{\pi }{6}n\right)}{5\pi n}\:=\:\frac{1}{30}\operatorname{sinc}\left(\frac{n}{6}\right)
\end{equation}
Therefore, the Fourier Transform of each would be
\begin{equation}
X_1\left(e^{j\omega }\right)\:=\:1\:\:\:\:0<\omega <\frac{\pi }{4}
\end{equation}
\begin{equation}
X_2\left(e^{j\omega }\right)\:=\:1\:\:\:\:0<\omega <\frac{\pi }{6}
\end{equation}
so using Parseval's theorem
\begin{equation}
S\:=\frac{1}{2\pi }\int _0^{2\pi }X_1\left(e^{j\omega }\right)X_2\left(e^{j\omega }\right)d\omega \:=\:\frac{1}{2\pi }\int _0^{\frac{\pi }{6}\:}\left(1\right)\left(1\right)d\omega \:\:=\:\frac{1}{12}
\end{equation}
Is this correct? am kinda doubting I did a mistake in the theory somewhere. Thanks in advance.


texts: Fourier Series Representation of Continuous-Time Periodic Signals

texts:
As a novice in signal processing, I have been going through Signals & Systems by Oppenheim to try and understand how continuous time periodic signals are represented by Fourier series coefficients.

The book defines a pair of equations which represent the Fourier series of a periodic continuous time signal.

I cannot understand when one should use the synthesis equation (3.38) and when one should use the analysis equation (3.39) to  determine the set of Fourier series coefficients.

Can someone give examples/explain as to when I should use one equation over another? Why is this so?


texts: If I want to apply a filter to an Image should I use the IIR or FIR filters?

texts: So I have discrete pixel values in the image and each row or column would make a signal. I would like to filter some of the rows or columns to remove high frequency components from them. Should I use the Z transform (difference equations) or work with the FFT of the data to perform the filtering?


texts: system response: time vs frequency. Why do I get different magnitudes?

texts:
I am simulating a linear system described by the following equation

$m \ddot{x} + c \dot{x} + k x = A \sin(\omega t)$

Fs = 1000;      % sampling frequency
tspan = 0:1/Fs:10;
m = 25;
c = 15;
k = 330000;
A = 100;           % forcing amplitude
omega = 114.89;    % forcing frequency

% State space function
odefun = @(t,u)[u(2); 1/m*(-k*u(1)-c*u(2)+A*sin(omega*t))];

[time,u] = ode45(odefun,tspan,[0 0]);

spost = u(:,1);    % displacement
vel = u(:,2);      % velocity

Then I compute the FFT of the displacement:

L = length(spost);
NFFT = 2^nextpow2(L);
SPOST = fft(spost,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);

And this is what I get:

As you can see the maximum value of the displacement is different in the 2 cases, but this does not look right to me. In fact, if I run a structural analysis in the time or in the frequency domain I would expect to obtain the same maximum displacement (since the excitation frequency is equal to the resonance frequency).

How can this be explained?

EDIT

I forgot to say I tried also increasing the simulation time, but still the 2 values do not match